 /*	THIS ESTIMATES DEMAND FOR OILS USING HAUSMAN INSTRUMENTS.  IN THIS VARIANT, EACH PRICE IS REGRESSED ON ALL OTHER REGIONS PRICES IN THE FIRST STAGE.  THERE COEFFICIENTS ON EACH REGION'S PRICE IS UNCONSTRAINED*/
set more off
#delimit;
cap log close;
*log using check, replace;
clear all;
set matsize 1000;
set mem 1200m;
set seed 10161978;


local path "c:/data/restat_12_9/analysis_data";
local estimator "newey2";
local ivestimator "newey2";
local lags "lag(4)";
local opt "se 3aster bd(2) td(2)";

local nbrands1=5;
local nbrands2=4;
local nbrands=9;
local nregions=47;
local list1 "AUNT_JEMIMA HUNGRY_JACK LOG_CABIN MRS_BUTTERWORTH PRIVATE_LABEL";
local list2 "AUNT_JEMIMA_LIGHT JACK_LITE LOG_CABIN_LIGHT MRS_B_LITE";


di "`path'";
do "c:/data/restat_12_9/programs/csolve_do/csolve";
do "c:/data/restat_12_9/programs/csolve_do/callf";
do "c:/data/restat_12_9/programs/csolve_do/callf_setup";
do "c:/data/restat_12_9/programs/csolve_do/syraids2";
use "`path'/syrmaster";



/*generate a time variable*/
sort REGION BRAND year month day;
bysort REGION BRAND: gen T=_n;

/*id variable for each brand/region pair*/
sort T REGION BRAND;
bysort T: gen id=_n;

/*make month and year dummies*/
qui tab month, gen(M);
qui tab year, gen(year);


/*make list of log price interacted with brand dummies for each segment. this will be used in AIDS.  lpreg1_`x'_`y' is price of brand `y' interacted with dummy variable if observation is brand`x'.
exclude interaction with pl because of adding up.  same for instruments*/
foreach y of local list1{;
	forval x=1/4{;
		local reglist1 "`reglist1' lpreg_`x'_`y'";
		local ivreglist1 "`ivreglist1' ivlpreg_`x'_`y'";
	};
};

foreach y of local list2{;
	forval x=1/3{;
		local reglist2 "`reglist2' lplight_`x'_`y'";
		local ivreglist2 "`ivreglist2' ivlplight_`x'_`y'";
	};
};


/*this creates 46 region dummies interacted with brand dummies*/
forval x=2/47{;
	forval y=1/5{;
		gen fe1`x'_`y'=reg`x'*bra1`y';
	};
};
forval x=2/47{;
	forval y=1/4{;
		gen fe2`x'_`y'=reg`x'*bra2`y';
	};
};


/*list of 46 region dummies, one for each of the share equations*/
forval x=2/47{;
	forval y=1/4{;
		local regionxbrand1 "`regionxbrand1' fe1`x'_`y'";
	};
};

forval x=2/47{;
	forval y=1/3{;
		local regionxbrand2 "`regionxbrand2' fe2`x'_`y'";
	};
};

forval x=1/8{;
	forval y=1/5{;		
		gen M1`x'b`y'=M`x'*bra1`y';
	};
};

forval x=1/8{;
	forval y=1/4{;		
		gen M2`x'b`y'=M`x'*bra2`y';
	};
};

forval x=1/4{;
	gen ivlaidpreg`x'=ivlaidpreg*bra1`x';
};

forval x=1/3{;
	gen ivlaidplight`x'=ivlaidplight*bra2`x';
};



/*this creates month dummies for each share equation*/
forval x=1/8{;
	forval y=1/4{;
		local Mlist1 "`Mlist1' M1`x'b`y'";
	};
};

forval x=1/8{;
	forval y=1/3{;
		local Mlist2 "`Mlist2' M2`x'b`y'";
	};
};


/*this creates time trends for each of 4 share equations*/
forval x=1/5{;
	gen t1`x'=T*bra1`x';
};

forval x=1/4{;
	gen t2`x'=T*bra2`x';
};





/*UNCONSTRAINED MODEL*/
encode REGION, gen(reg);
tsset id T;
qui compress;

/*AIDS OLS ESTIMATES FOR THE REGULAR SEGMENT*/
`estimator' share laidpreg1 laidpreg2 laidpreg3 laidpreg4 `reglist1' `Mlist1' `regionxbrand1' bra11-bra14 t11-t14 if BRAND~="PRIVATE LABEL" & nest==1, nocon `lags';
 
/*store point estimate and vc matrix*/
matrix aidB=e(b);
matrix aidVC=e(V);

/*coefficients on log quantity term*/
forval x=1/4{;
	scalar b`x'=_coef[laidpreg`x'];
};

/*recover Private Label share equation parameters through adding up restrictions of AIDS model*/
scalar b5=-b1-b2-b3-b4;

/*g`x'p`y' is the coefficient on the (log) price of brand `y' in equation `x'*/
forval x=1/4{;
	foreach y of local list1{;
		scalar g`x'p`y'=_coef[lpreg_`x'_`y'];
	};
};
foreach y of local list1{;
	scalar g5p`y'=-_coef[lpreg_1_`y']-_coef[lpreg_2_`y']-_coef[lpreg_3_`y']-_coef[lpreg_4_`y'];
};
/*rename g's*/
forval x=1/5{;
	scalar g`x'p1=g`x'pAUNT_JEMIMA;
	scalar g`x'p2=g`x'pHUNGRY_JACK;
	scalar g`x'p3=g`x'pLOG_CABIN;
	scalar g`x'p4=g`x'pMRS_BUTTERWORTH;
	scalar g`x'p5=g`x'pPRIVATE_LABEL;
};




/*Average share by brand in regular segment, to be used in evaluating elasticities at grand means.*/
forval x=1/`nbrands1'{;
	sum share if bra1`x'==1;
	scalar s1`x'=r(mean);
};

/*Test for stability across regions*/
preserve;
	gen west=0;
	replace west=1 if REGION=="DALLS"|REGION=="DENVER"|REGION=="DESMOINES"|REGION=="HOUSTON"|REGION=="KC"|REGION=="LROCK"|REGION=="MINN"|REGION=="OKCITY"|REGION=="OMAHA"|REGION=="PHOENIX"|REGION=="SANANT"|REGION=="STL"|REGION=="WTEX";
	gen laidpreg_w=laidpreg*west;
	gen lpriceAUNT_JEMIMA_w=lpriceAUNT_JEMIMA*west;
	gen lpriceHUNGRY_JACK_w= lpriceHUNGRY_JACK*west;
	gen  lpriceLOG_CABIN_w= lpriceLOG_CABIN*west;
	gen  lpriceMRS_BUTTERWORTH_w= lpriceMRS_BUTTERWORTH*west;
	gen  lpricePRIVATE_LABEL_w= lpricePRIVATE_LABEL*west;
	
	newey2 share west laidpreg laidpreg_w lpriceAUNT_JEMIMA lpriceAUNT_JEMIMA_w lpriceHUNGRY_JACK lpriceHUNGRY_JACK_w lpriceLOG_CABIN lpriceLOG_CABIN_w lpriceMRS_BUTTERWORTH lpriceMRS_BUTTERWORTH_w lpricePRIVATE_LABEL lpricePRIVATE_LABEL_w M2-M9 if BRAND=="LOG_CABIN", lag(4);
	test (laidpreg_w=0) (lpriceAUNT_JEMIMA_w=0) (lpriceHUNGRY_JACK_w=0) (lpriceLOG_CABIN_w=0) (lpriceMRS_BUTTERWORTH_w=0) (lpricePRIVATE_LABEL_w=0) (west=0);
	newey2 share west laidpreg laidpreg_w lpriceAUNT_JEMIMA lpriceAUNT_JEMIMA_w lpriceHUNGRY_JACK lpriceHUNGRY_JACK_w lpriceLOG_CABIN lpriceLOG_CABIN_w lpriceMRS_BUTTERWORTH lpriceMRS_BUTTERWORTH_w lpricePRIVATE_LABEL lpricePRIVATE_LABEL_w M2-M9 if BRAND=="MRS_BUTTERWORTH", lag(4);
	test (laidpreg_w=0) (lpriceAUNT_JEMIMA_w=0) (lpriceHUNGRY_JACK_w=0) (lpriceLOG_CABIN_w=0) (lpriceMRS_BUTTERWORTH_w=0) (lpricePRIVATE_LABEL_w=0) (west=0);
restore;
	/*IV AIDS ESTIMATES FOR THE REGULAR SEGMENT*/
/*First stage diagnostics.  ESTIMATE WITH SINGLE EQUATIONS AND COMPARE TO MAKE SURE THIS IS RIGHT.*/


*ivregress 2sls share `Mlist1' `regionxbrand1' bra11-bra14 count11-count14  (laidpreg1 laidpreg2 laidpreg3 laidpreg4 `reglist1'=ivlaidpreg1 ivlaidpreg2 ivlaidpreg3 ivlaidpreg4 `ivreglist1') if BRAND~="PRIVATE LABEL" & nest==1, nocon `lags';
*estat firststage;

`ivestimator' share `Mlist1' `regionxbrand1' bra11-bra14 t11-t14  (laidpreg1 laidpreg2 laidpreg3 laidpreg4 `reglist1'=ivlaidpreg1 ivlaidpreg2 ivlaidpreg3 ivlaidpreg4 `ivreglist1') if BRAND~="PRIVATE LABEL" & nest==1, nocon `lags';



matrix ivaidB=e(b);
matrix ivaidVC=e(V);
forval x=1/4{;
	scalar ivb`x'=_coef[laidpreg`x'];
};
/*recover Private Label share equation parameters through adding up restrictions of AIDS model*/
scalar ivb5=-ivb1-ivb2-ivb3-ivb4;
#delimit;
forval x=1/4{;
	foreach y of local list1{;
		scalar ivg`x'p`y'=_coef[lpreg_`x'_`y'];
	};
};
foreach y of local list1{;
	scalar ivg5p`y'=-_coef[lpreg_1_`y']-_coef[lpreg_2_`y']-_coef[lpreg_3_`y']-_coef[lpreg_4_`y'];
};

forval x=1/5{;
	scalar ivg`x'p1=ivg`x'pAUNT_JEMIMA;
	scalar ivg`x'p2=ivg`x'pHUNGRY_JACK;
	scalar ivg`x'p3=ivg`x'pLOG_CABIN;
	scalar ivg`x'p4=ivg`x'pMRS_BUTTERWORTH;
	scalar ivg`x'p5=ivg`x'pPRIVATE_LABEL;
};




/*OLS AIDS ESTIMATES FOR THE LIGHT SEGMENT*/
`estimator' share laidplight1 laidplight2 laidplight3 `reglist2' `Mlist2' `regionxbrand2' bra21-bra23 t21-t23 if BRAND~="MRS_B_LITE" & nest==2, nocon `lags';
matrix aidB2=e(b);
matrix aidVC2=e(V);
forval x=1/3{;
	scalar b2`x'=_coef[laidplight`x'];
};

/*recover Private Label share equation parameters through adding up restrictions of AIDS model*/
scalar b24=-b21-b22-b23;
#delimit;
forval x=1/3{;
	foreach y of local list2{;
		scalar g2`x'p`y'=_coef[lplight_`x'_`y'];
	};
};
foreach y of local list2{;
	scalar g24p`y'=-_coef[lplight_1_`y']-_coef[lplight_2_`y']-_coef[lplight_3_`y'];
};

forval x=1/4{;
	scalar g2`x'p1=g2`x'pAUNT_JEMIMA_LIGHT;
	scalar g2`x'p2=g2`x'pJACK_LITE;
	scalar g2`x'p3=g2`x'pLOG_CABIN_LIGHT;
	scalar g2`x'p4=g2`x'pMRS_B_LITE;
};

forval x=1/`nbrands2'{;
	sum share if bra2`x'==1;
	scalar s2`x'=r(mean);
};

/*HERE ARE THE AIDS IV ESTIMATES FOR THE LIGHT SEGMENT*/
*ivregress 2sls share `Mlist2' `regionxbrand2' bra21-bra23 count21-count23 (laidplight1 laidplight2 laidplight3 `reglist2'=ivlaidplight1 ivlaidplight2 ivlaidplight3`ivreglist2' ) /*reg1-reg47*/  if BRAND~="MRS_B_LITE" & nest==2, nocon `lags';
*estat firststage;


`ivestimator' share `Mlist2' `regionxbrand2' bra21-bra23 t21-t23 (laidplight1 laidplight2 laidplight3 `reglist2'=ivlaidplight1 ivlaidplight2 ivlaidplight3`ivreglist2' ) /*reg1-reg47*/  if BRAND~="MRS_B_LITE" & nest==2, nocon `lags';

matrix ivaidB2=e(b);
matrix ivaidVC2=e(V);
forval x=1/3{;
	scalar ivb2`x'=_coef[laidplight`x'];
};

/*recover Private Label share equation parameters through adding up restrictions of AIDS model*/
scalar ivb24=-ivb21-ivb22-ivb23;
#delimit;
forval x=1/3{;
	foreach y of local list2{;
		scalar ivg2`x'p`y'=_coef[lplight_`x'_`y'];
	};
};
foreach y of local list2{;
	scalar ivg24p`y'=-_coef[lplight_1_`y']-_coef[lplight_2_`y']-_coef[lplight_3_`y'];
};

forval x=1/4{;
	scalar ivg2`x'p1=ivg2`x'pAUNT_JEMIMA_LIGHT;
	scalar ivg2`x'p2=ivg2`x'pJACK_LITE;
	scalar ivg2`x'p3=ivg2`x'pLOG_CABIN_LIGHT;
	scalar ivg2`x'p4=ivg2`x'pMRS_B_LITE;
};



/*Hosken's idea: weight volume and expenditures by population in linear model.*/
gen wvol=vol/POP;
forval x=1/5{;
	gen wexpreg`x'=expreg`x'/POP;
};

/*Create list of average prices interacted with brand dummies for linear equations*/
foreach y of local list1{;
	forval x=1/5{;
		local linplistreg "`linplistreg' preg_`x'_`y'";
	};
};
foreach y of local list1{;
	forval x=1/5{;
		local ivlinplistreg "`ivlinplistreg' ivpreg_`x'_`y'";
	};
};

/*List of region fixed effects for each brand's equation*/
forval x=2/47{;
	forval y=1/5{;
		local lregionxbrand "`lregionxbrand' fe1`x'_`y'";
	};
};

/*Month dummies interacted with brand dummies*/
local lMlist1 "`Mlist1' M11b5 M12b5 M13b5 M14b5 M15b5 M16b5 M17b5 M18b5";

/*OLS estimates of linear demand, regular segment*/
`estimator' wvol t11-t15 bra11-bra15 `lMlist1' `lregionxbrand' wexpreg1-wexpreg5 `linplistreg' if nest==1, nocon `lags';  

matrix linB=e(b);
matrix linVC=e(V);
#delimit;
/*p_i_j is coefficient on price j in equation i*/
forval x=1/5{;
	scalar p_`x'_1=_coef[preg_`x'_AUNT_JEMIMA];
	scalar p_`x'_2=_coef[preg_`x'_HUNGRY_JACK];
	scalar p_`x'_3=_coef[preg_`x'_LOG_CABIN];
	scalar p_`x'_4=_coef[preg_`x'_MRS_BUTTERWORTH];
	scalar p_`x'_5=_coef[preg_`x'_PRIVATE_LABEL];
};


forval x=1/5{;
    forval y=1/5{;
        scalar linb`x'p`y'=p_`x'_`y';
    };
};

mat B_reg=J(5,5,0);
forval x=1/5{;
    forval y=1/5{;
        mat B_reg[`x',`y']=linb`x'p`y';
    };
};

mat theta_reg=J(5,1,0);

forval x=1/5{;
	mat theta_reg[`x',1]=_coef[wexpreg`x'];
};	


/*Average prices and wvol for use in evaluating elasticities at grand means*/
forval x=1/5{;
	sum price if bra1`x'==1;
	scalar p1`x'=r(mean);
};

forval x=1/`nbrands1'{;
	sum wvol if bra1`x'==1;
	scalar w1`x'=r(mean);
};

/*2sls estimates of linear equation*/
`ivestimator' wvol t11-t15 bra11-bra15 `lMlist1' `lregionxbrand' (wexpreg1-wexpreg5 `linplistreg'=ivwexpreg1-ivwexpreg5 `ivlinplistreg') if nest==1, nocon `lags';  


matrix ivlinB=e(b);
matrix ivlinVC=e(V);
#delimit;
forval x=1/5{;
	scalar ivp_`x'_1=_coef[preg_`x'_AUNT_JEMIMA];
	scalar ivp_`x'_2=_coef[preg_`x'_HUNGRY_JACK];
	scalar ivp_`x'_3=_coef[preg_`x'_LOG_CABIN];
	scalar ivp_`x'_4=_coef[preg_`x'_MRS_BUTTERWORTH];
	scalar ivp_`x'_5=_coef[preg_`x'_PRIVATE_LABEL];
};


forval x=1/5{;
    forval y=1/5{;
        scalar ivlinb`x'p`y'=ivp_`x'_`y';
    };
};

mat ivB_reg=J(5,5,0);
forval x=1/5{;
    forval y=1/5{;
        mat ivB_reg[`x',`y']=ivlinb`x'p`y';
    };
};


mat ivtheta_reg=J(5,1,0);

forval x=1/5{;
	mat ivtheta_reg[`x',1]=_coef[wexpreg`x'];
};	


/*population weighted expenditure term for light syrups*/
forval x=1/4{;
	gen wexplight`x'=explight`x'/POP;
};

/*lists of prices for light syrups*/
foreach y of local list2{;
	forval x=1/4{;
		local linplistlight "`linplistlight' plight_`x'_`y'";
	};
};

foreach y of local list2{;
	forval x=1/4{;
		local ivlinplistlight "`ivlinplistlight' ivplight_`x'_`y'";
	};
};


forval x=2/47{;
	forval y=1/4{;
		local lregionxbrand2 "`lregionxbrand2' fe2`x'_`y'";
	};
};

local lMlist2 "`Mlist2' M21b4 M22b4 M23b4 M24b4 M25b4 M26b4 M27b4 M28b4";


/*ols estimates of linear demand for light segment*/
`estimator' wvol t21-t24 bra21-bra24 `lMlist2' `lregionxbrand2' wexplight1-wexplight4 `linplistlight' if nest==2, nocon `lags';  

matrix linB2=e(b);
matrix linVC2=e(V);
#delimit;
forval x=1/4{;
	scalar p2_`x'_1=_coef[plight_`x'_AUNT_JEMIMA_LIGHT];
	scalar p2_`x'_2=_coef[plight_`x'_JACK_LITE];
	scalar p2_`x'_3=_coef[plight_`x'_LOG_CABIN_LIGHT];
	scalar p2_`x'_4=_coef[plight_`x'_MRS_B_LITE];
};


forval x=1/4{;
    forval y=1/4{;
        scalar linb2`x'p`y'=p2_`x'_`y';
    };
};

mat B_light=J(4,4,0);
forval x=1/4{;
    forval y=1/4{;
        mat B_light[`x',`y']=linb2`x'p`y';
    };
};

mat theta_light=J(4,1,0);

forval x=1/4{;
	mat theta_light[`x',1]=_coef[wexplight`x'];
};	

/*average prices and wvol for evaluating elasticities at grand means*/
forval x=1/4{;
	sum price if bra2`x'==1;
	scalar p2`x'=r(mean);
};

forval x=1/`nbrands2'{;
	sum wvol if bra2`x'==1;
	scalar w2`x'=r(mean);
};

/*2sls estimates of linear demand for light segment*/
`ivestimator' wvol t21-t24 bra21-bra24 `lMlist2' `lregionxbrand2'  (wexplight1-wexplight4 `linplistlight'=ivwexplight1-ivwexplight4 `ivlinplistlight') if nest==2, nocon `lags';  


matrix ivlinB2=e(b);
matrix ivlinVC2=e(V);
#delimit;
forval x=1/4{;
	scalar ivp2_`x'_1=_coef[plight_`x'_AUNT_JEMIMA_LIGHT];
	scalar ivp2_`x'_2=_coef[plight_`x'_JACK_LITE];
	scalar ivp2_`x'_3=_coef[plight_`x'_LOG_CABIN_LIGHT];
	scalar ivp2_`x'_4=_coef[plight_`x'_MRS_B_LITE];
};


forval x=1/4{;
    forval y=1/4{;
        scalar ivlinb2`x'p`y'=ivp2_`x'_`y';
    };
};

mat ivB_light=J(4,4,0);
forval x=1/4{;
    forval y=1/4{;
        mat ivB_light[`x',`y']=ivlinb2`x'p`y';
    };
};

mat ivtheta_light=J(4,1,0);

forval x=1/4{;
	mat ivtheta_light[`x',1]=_coef[wexplight`x'];
};	



save temp, replace;


/*Create a dataset where an observation is a region/time/segment market.  Use this to estimate middle level.*/  
drop vol-APRICE merged-lprice lexpreg-ivprice bra11-reg47 ivlaidpreg-ivlaidplight laidpreg1-ivplight_4_MRS_B_LITE ivwexpreg-ivwexplight ivwexpreg1-ivwexplight4 BRAND livprice id ivlaidpreg1-ivlaidplight3 wvol wexpreg1-wexplight4 t11-t24 fe12_1-M28b4;
duplicates drop;


/*Middle segment expenditures.  This is (regular expenditure)/POP+(light expenditure)/POP*/
gen wexpmid=wexpreg+wexplight;

/*"real" expenditures.  This will be dependent variable in middle level.*/
gen lnvolmid=0;
replace lnvolmid=laidpreg if nest==1;
replace lnvolmid=laidplight if nest==2;
 
/*gen nest dummies*/
qui tab nest, gen(nest);
/*Nest dummies interacted with log price index.  lpx_y is (log) price index of segment x in middle segment y.*/
gen lp1_1=nest1*indexreg;
gen lp1_2=nest1*indexlight;
gen lp2_1=nest2*indexreg; 
gen lp2_2=nest2*indexlight;

gen ivlp1_1=nest1*ivindexreg;
gen ivlp1_2=nest1*ivindexlight;
gen ivlp2_1=nest2*ivindexreg;
gen ivlp2_2=nest2*ivindexlight;

/*Nest dummies interacted with price index for linear middle level.  px_y is (log) price index of segment x in middle segment y.*/
gen p1_1=nest1*lin_indexreg;
gen p1_2=nest1*lin_indexlight;
gen p2_1=nest2*lin_indexreg;
gen p2_2=nest2*lin_indexlight;

gen ivp1_1=nest1*ivlin_indexreg;
gen ivp1_2=nest1*ivlin_indexlight;
gen ivp2_1=nest2*ivlin_indexreg;
gen ivp2_2=nest2*ivlin_indexlight;


/*Log expenditures interacted with nest dummies*/
gen lexp1=lexp*nest1;
gen lexp2=lexp*nest2;


/* expenditures/population interacted with nest dummies*/
gen wexp1=wexp*nest1;
gen wexp2=wexp*nest2;

gen ivlexp1=ivlexp*nest1;
gen ivlexp2=ivlexp*nest2;
gen ivwexp1=ivwexp*nest1;
gen ivwexp2=ivwexp*nest2;


qui tab REGION, gen(reg);

/*Nest*region dummies*/
forval x=2/47{;
	forval y=1/2{;
		gen fe`y'_`x'=nest`y'*reg`x';
	};
};

/*Nest*month dummies*/
forval x=1/8{;
	forval y=1/2{;		
		gen M`x'b`y'=M`x'*nest`y';
	};
};

#delimit;
forval x=2/47{;
	forval y=1/2{;
		local femid "`femid' fe`y'_`x'";
	};
};

/*this creates month dummies for each of 2 middle level equations*/
forval x=1/8{;
	forval y=1/2{;
		local Mlistmid "`Mlistmid' M`x'b`y'";
	};
};

gen count1=t*nest1;
gen count2=t*nest2;
sort nest REGION t;


sort t REGION nest;
/*separate id for each nest/region combo*/
bysort t: gen id=_n;
tsset id T;

/*Middle Level for Aids Model*/
`estimator' lnvolmid `femid' `Mlistmid' nest1 nest2 lexp1 lexp2 lp1_1 lp1_2 lp2_1 lp2_2 count1 count2, nocon `lags';
matrix aidMidB=e(b);
matrix aidMidVC=e(V);


sca delta_mid_reg_reg=_coef[lp1_1];
sca delta_mid_reg_light=_coef[lp1_2];
sca delta_mid_light_reg=_coef[lp2_1];
sca delta_mid_light_light=_coef[lp2_2];
sca beta_reg=_coef[lexp1];
sca beta_light=_coef[lexp2];

mat deltamid=[delta_mid_reg_reg, delta_mid_reg_light\delta_mid_light_reg, delta_mid_light_light];
mat betamid=[beta_reg\beta_light];


`ivestimator' lnvolmid `femid' `Mlistmid' nest1 nest2 count1 count2 (lexp1 lexp2 lp1_1 lp1_2 lp2_1 lp2_2=ivlexp1 ivlexp2 ivlp1_1 ivlp1_2 ivlp2_1 ivlp2_2), nocon `lags';
sca ivdelta_mid_reg_reg=_coef[lp1_1];
sca ivdelta_mid_reg_light=_coef[lp1_2];
sca ivdelta_mid_light_reg=_coef[lp2_1];
sca ivdelta_mid_light_light=_coef[lp2_2];
sca ivbeta_reg=_coef[lexp1];
sca ivbeta_light=_coef[lexp2];

mat ivdeltamid=[ivdelta_mid_reg_reg, ivdelta_mid_reg_light\ivdelta_mid_light_reg, ivdelta_mid_light_light];
mat ivbetamid=[ivbeta_reg\ivbeta_light];

/*For linear model*/
`estimator' wexpmid `femid' `Mlistmid' nest1 nest2 wexp1 wexp2 count1 count2 p1_1 p1_2 p2_1 p2_2, nocon `lags';
sca ldelta_mid_reg_reg=_coef[p1_1];
sca ldelta_mid_reg_light=_coef[p1_2];
sca ldelta_mid_light_reg=_coef[p2_1];
sca ldelta_mid_light_light=_coef[p2_2];
sca lbeta_reg=_coef[wexp1];
sca lbeta_light=_coef[wexp2];

mat ldeltamid=[ldelta_mid_reg_reg, ldelta_mid_reg_light\ldelta_mid_light_reg, ldelta_mid_light_light];
mat lbetamid=[lbeta_reg\lbeta_light];


`ivestimator' wexpmid `femid' `Mlistmid' nest1 nest2 count1 count2 (wexp1 wexp2 p1_1 p1_2 p2_1 p2_2=ivwexp1 ivwexp2 ivp1_1 ivp1_2 ivp2_1 ivp2_2), nocon `lags';
sca livdelta_mid_reg_reg=_coef[p1_1];
sca livdelta_mid_reg_light=_coef[p1_2];
sca livdelta_mid_light_reg=_coef[p2_1];
sca livdelta_mid_light_light=_coef[p2_2];
sca livbeta_reg=_coef[wexp1];
sca livbeta_light=_coef[wexp2];

mat livdeltamid=[livdelta_mid_reg_reg, livdelta_mid_reg_light\livdelta_mid_light_reg, livdelta_mid_light_light];
mat livbetamid=[livbeta_reg\livbeta_light];


drop nest lexp-laidplight lin_indexreg lin_indexlight ivlexp ivindexreg ivindexlight ivlin_indexreg ivlin_indexlight ivexp wexpreg-ivwexp2 fe1_2- M8b2 count1 count2 id;

duplicates drop;

gen i=0;
forval x=1/47{;
	replace i=`x' if reg`x'==1;
};

sort i t;

bysort i: gen count=_n;
tsset i count;
qui tab month, gen(M);
/*2sls Top For Aids Model*/
`ivestimator' lqall  M2-M9 reg2-reg47 count (indexall=ivindexall), `lags';
scalar ivdelta_top=_coef[indexall];
matrix iv_b_top=e(b);
matrix iv_vc_top=e(V);
mata: iv_b_top=st_matrix("iv_b_top");
mata: iv_vc_top=st_matrix("iv_vc_top");



/*ols Top For Aids Model*/
`estimator' lqall indexall M2-M9 count reg2-reg47, `lags';
scalar delta_top=_coef[indexall];
matrix ols_b_top=e(b);
matrix ols_vc_top=e(V);
mata: ols_b_top=st_matrix("ols_b_top");
mata: ols_vc_top=st_matrix("ols_vc_top");


mata: mata matsave "c:/data/restat_12_9/analysis_data/top_level" ols_b_top ols_vc_top iv_b_top iv_vc_top, replace;

/*For linear model*/
`ivestimator' wexp  M2-M9 reg2-reg47 count (lin_indexall=ivlin_indexall), `lags';
mat ivldelta_top=_coef[lin_indexall];
`estimator' wexp lin_indexall M2-M9 count reg2-reg47, `lags';
mat ldelta_top=_coef[lin_indexall];



#delimit;
use temp, replace;
gen bra1=0;
replace bra1=1 if BRAND=="AUNT_JEMIMA";
gen bra2=0;
replace bra2=1 if BRAND=="HUNGRY_JACK";
gen bra3=0;
replace bra3=1 if BRAND=="LOG_CABIN";
gen bra4=0;
replace bra4=1 if BRAND=="MRS_BUTTERWORTH";
gen bra5=0;
replace bra5=1 if BRAND=="PRIVATE_LABEL";
gen bra6=0;
replace bra6=1 if BRAND=="AUNT_JEMIMA_LIGHT";
gen bra7=0;
replace bra7=1 if BRAND=="JACK_LITE";
gen bra8=0;
replace bra8=1 if BRAND=="LOG_CABIN_LIGHT";
gen bra9=0;
replace bra9=1 if BRAND=="MRS_B_LITE";




gen expmid=expreg+explight;


#delimit;
/*total syrup sales by type over time*/
sort BRAND REGION nest t;
bysort BRAND REGION nest : egen tsaletype=total(expmid);

/*total syrup sales over time*/
bysort BRAND REGION: egen tsaleall=total(exp);

/*total syrup sales by brand over time*/
bysort BRAND REGION: egen bsale=total(sale);
gen avgsharenest=bsale/tsaletype;
gen avgshare=bsale/tsaleall;


 /*w`x'b`y'reg is the average revenue share of producy `y' in region `x' conditional on regular*/
forval x=1/`nregions'{;
    forval y=1/`nbrands1'{;	    
        qui sum avgsharenest if reg`x'==1 & bra1`y'==1;
        scalar w`x'b`y'reg=r(mean);
    };
};

/*w`x'b`y'light is the average revenue share of producy `y' in region `x' conditional on light*/
forval x=1/`nregions'{;
    forval y=1/`nbrands2'{;
        qui sum avgsharenest if reg`x'==1 & bra2`y'==1;
        scalar w`x'b`y'light=r(mean);
    };
};

/*w`x'b`y'all is the average revenue share of producy `y' in region `x'*/
forval x=1/`nregions'{;
    forval y=1/`nbrands'{;
        qui sum avgshare if reg`x'==1 & brand`y'==1;
        scalar w`x'b`y'all=r(mean);
    };
};




/*total sales by segment over all regions and time*/
bysort BRAND nest: egen national_tsale_type=total(expmid);

/*total sales over all regions and time*/
bysort BRAND: egen national_tsale_all=total(exp);


bysort BRAND: egen national_bsale=total(sale);
gen nat_avgshare_nest=national_bsale/national_tsale_type;
gen nat_avgshare_all=national_bsale/national_tsale_all;




/*national sales shares by brand within nest*/
forval x=1/`nbrands1'{;
	sum nat_avgshare_nest if bra1`x'==1;
	scalar avgw`x'reg=r(mean);
};

forval x=1/`nbrands2'{;
	sum nat_avgshare_nest if bra2`x'==1;
	scalar avgw`x'light=r(mean);
};



/*uncondtional national sales shares by brand*/
forval x=1/`nbrands1'{;
	sum nat_avgshare_all if bra1`x'==1;
	scalar avgw`x'reg_all=r(mean);
};

forval x=1/`nbrands2'{;
	sum nat_avgshare_all if bra2`x'==1;
	scalar avgw`x'light_all=r(mean);
};



/*HERE WE CONSTRUCT OVER-ALL ELASTICITIES EVALUATED AT AVERAGES TAKEN OVER REGIONS AND TIME*/


/*AIDS estimated by OLS*/
/*regular*/
/*within same nest elasticities*/
forval x=1/`nbrands1'{;
	forval y=1/`nbrands1'{;
		scalar ereg`x'`y'=(1+delta_mid_reg_reg)*avgw`y'reg+beta_reg*(1+delta_top)*avgw`y'reg_all+g`x'p`y'/s1`x'+(b`x'/s1`x')*(delta_mid_reg_reg*avgw`y'reg+
		avgw`y'reg_all*(1+delta_top)*beta_reg);
	};
};
forval x=1/`nbrands2'{;
	forval y=1/`nbrands2'{;
		scalar elight`x'`y'=(1+delta_mid_light_light)*avgw`y'light+beta_light*(1+delta_top)*avgw`y'light_all+g2`x'p`y'/s2`x'+(b2`x'/s2`x')*(delta_mid_light_light*avgw`y'light+
		avgw`y'light_all*(1+delta_top)*beta_light);
	};
};


forval x=1/`nbrands1'{;
	scalar ereg`x'`x'=-1+ereg`x'`x';
};
forval x=1/`nbrands2'{;
	scalar elight`x'`x'=-1+elight`x'`x';
};
mat regnestelas=J(5,5,0);
mat lightnestelas=J(4,4,0);


forval x=1/5{;
	forval y=1/5{;
	mat regnestelas[`x',`y']=ereg`x'`y';
	};
};

forval x=1/4{;
	forval y=1/4{;
	mat lightnestelas[`x',`y']=elight`x'`y';
	};
};

/*Cross nest elasticities*/
forval x=1/5{;
	forval y=1/4{;
		scalar ereg`x'light`y'=(1+b`x'/s1`x')*(beta_reg*(1+delta_top)*avgw`y'light_all+delta_mid_reg_light*avgw`y'light);
	};
};


mat reg_lightelas=J(5,4,0);

forval x=1/5{;
	forval y=1/4{;
	mat reg_lightelas[`x',`y']=ereg`x'light`y';
	};
};
	

forval x=1/4{;
	forval y=1/5{;
		scalar elight`x'reg`y'=(1+b2`x'/s2`x')*(beta_light*(1+delta_top)*avgw`y'reg_all+delta_mid_light_reg*avgw`y'reg);
	};
};
mat light_regelas=J(4,5,0);

forval x=1/4{;
	forval y=1/5{;
	mat light_regelas[`x',`y']=elight`x'reg`y';
	};
};


/*Elasticities evaluated at overall means for syrup AIDS model estimated with OLS*/
mat aidselasticities=J(9,9,0);
mat aidselasticities[1,1]=regnestelas;
mat aidselasticities[6,6]=lightnestelas;
mat aidselasticities[1,6]=reg_lightelas;
mat aidselasticities[6,1]=light_regelas;


/*Within nest elasticities syrup AIDS model estimated with 2SLS*/
forval x=1/`nbrands1'{;
	forval y=1/`nbrands1'{;
		scalar ivereg`x'`y'=(1+ivdelta_mid_reg_reg)*avgw`y'reg+ivbeta_reg*(1+ivdelta_top)*avgw`y'reg_all+ivg`x'p`y'/s1`x'+(ivb`x'/s1`x')*(ivdelta_mid_reg_reg*avgw`y'reg+
		avgw`y'reg_all*(1+ivdelta_top)*ivbeta_reg);
	};
};



forval x=1/`nbrands2'{;
	forval y=1/`nbrands2'{;
		scalar ivelight`x'`y'=(1+ivdelta_mid_light_light)*avgw`y'light+ivbeta_light*(1+ivdelta_top)*avgw`y'light_all+ivg2`x'p`y'/s2`x'+(ivb2`x'/s2`x')*(ivdelta_mid_light_light*avgw`y'light+
		avgw`y'light_all*(1+ivdelta_top)*ivbeta_light);
	};
};


forval x=1/`nbrands1'{;
	scalar ivereg`x'`x'=-1+ivereg`x'`x';
};
forval x=1/`nbrands2'{;
	scalar ivelight`x'`x'=-1+ivelight`x'`x';
};


mat ivregnestelas=J(5,5,0);
mat ivlightnestelas=J(4,4,0);


forval x=1/5{;
	forval y=1/5{;
	mat ivregnestelas[`x',`y']=ivereg`x'`y';
	};
};

forval x=1/4{;
	forval y=1/4{;
	mat ivlightnestelas[`x',`y']=ivelight`x'`y';
	};
};

/*Cross nest elasticities*/
forval x=1/5{;
	forval y=1/4{;
		scalar ivereg`x'light`y'=(1+ivb`x'/s1`x')*(ivbeta_reg*(1+ivdelta_top)*avgw`y'light_all+ivdelta_mid_reg_light*avgw`y'light);
	};
};

mat ivreg_lightelas=J(5,4,0);

forval x=1/5{;
	forval y=1/4{;
	mat ivreg_lightelas[`x',`y']=ivereg`x'light`y';
	};
};
	

forval x=1/4{;
	forval y=1/5{;
		scalar ivelight`x'reg`y'=(1+ivb2`x'/s2`x')*(ivbeta_light*(1+ivdelta_top)*avgw`y'reg_all+ivdelta_mid_light_reg*avgw`y'reg);
	};
};
mat ivlight_regelas=J(4,5,0);

forval x=1/4{;
	forval y=1/5{;
	mat ivlight_regelas[`x',`y']=ivelight`x'reg`y';
	};
};


/*Elasticities evaluated at overall means for syrup AIDS model estimated with 2SLS*/
mat ivaidselasticities=J(9,9,0);
mat ivaidselasticities[1,1]=ivregnestelas;
mat ivaidselasticities[6,6]=ivlightnestelas;
mat ivaidselasticities[1,6]=ivreg_lightelas;
mat ivaidselasticities[6,1]=ivlight_regelas;



/*Linear Elasticities*/
/*within nest elasticities*/

forval x=1/`nbrands1'{;
	forval y=1/`nbrands1'{;
		scalar linereg`x'`y'=(p1`y'/w1`x')*(linb`x'p`y'+theta_reg[`x',1]*(ldeltamid[1,1]*avgw`y'reg+avgw`y'reg_all*lbetamid[1,1]*ldelta_top[1,1]));
	};
};
mat lregelas=J(5,5,0);
forval x=1/`nbrands1'{;
	forval y=1/`nbrands1'{;
		mat lregelas[`x',`y']=linereg`x'`y';
	};
};
 
forval x=1/`nbrands2'{;
	forval y=1/`nbrands2'{;
		scalar linelight`x'`y'=(p2`y'/w2`x')*(linb2`x'p`y'+theta_light[`x',1]*(ldeltamid[2,2]*avgw`y'light+avgw`y'light_all*lbetamid[2,1]*ldelta_top[1,1]));
	};
};
mat llightelas=J(4,4,0);
forval x=1/`nbrands2'{;
	forval y=1/`nbrands2'{;
		mat llightelas[`x',`y']=linelight`x'`y';
	};
};
 
/*across nest elasticities*/
forval x=1/`nbrands1'{;
	forval y=1/`nbrands2'{;
		scalar linereg`x'light`y'=(p2`y'/w1`x')*(theta_reg[`x',1]*(ldeltamid[1,2]*avgw`y'light+lbetamid[1,1]*ldelta_top[1,1]*avgw`y'light_all ));
	};
};
mat lreg_lightelas=J(5,4,0);

forval x=1/`nbrands1'{;
	forval y=1/`nbrands2'{;
		mat lreg_lightelas[`x',`y']=linereg`x'light`y';
	};
};



forval x=1/`nbrands2'{;
	forval y=1/`nbrands1'{;
		scalar linelight`x'reg`y'=(p1`y'/w2`x')*(theta_light[`x',1]*(ldeltamid[2,1]*avgw`y'reg+lbetamid[2,1]*ldelta_top[1,1]*avgw`y'reg_all ));
	};
};

mat llight_regelas=J(4,5,0);

forval x=1/`nbrands2'{;
	forval y=1/`nbrands1'{;
		mat llight_regelas[`x',`y']=linelight`x'reg`y';
	};
};
/*Linear elasticities estimated by OLS*/
mat linelasticities=J(9,9,0);
mat linelasticities[1,1]=lregelas;
mat linelasticities[6,6]=llightelas;
mat linelasticities[1,6]=lreg_lightelas;
mat linelasticities[6,1]=llight_regelas;


forval x=1/`nbrands1'{;
	forval y=1/`nbrands1'{;
		scalar ivlinereg`x'`y'=(p1`y'/w1`x')*(ivlinb`x'p`y'+ivtheta_reg[`x',1]*(livdeltamid[1,1]*avgw`y'reg+avgw`y'reg_all*livbetamid[1,1]*ivldelta_top[1,1]));
	};
};
mat ivlregelas=J(5,5,0);
forval x=1/`nbrands1'{;
	forval y=1/`nbrands1'{;
		mat ivlregelas[`x',`y']=ivlinereg`x'`y';
	};
};
 
forval x=1/`nbrands2'{;
	forval y=1/`nbrands2'{;
		scalar ivlinelight`x'`y'=(p2`y'/w2`x')*(ivlinb2`x'p`y'+ivtheta_light[`x',1]*(livdeltamid[2,2]*avgw`y'light+avgw`y'light_all*livbetamid[2,1]*ivldelta_top[1,1]));
	};
};
mat ivllightelas=J(4,4,0);
forval x=1/`nbrands2'{;
	forval y=1/`nbrands2'{;
		mat ivllightelas[`x',`y']=ivlinelight`x'`y';
	};
};
 
/*across nest elasticities*/
forval x=1/`nbrands1'{;
	forval y=1/`nbrands2'{;
		scalar ivlinereg`x'light`y'=(p2`y'/w1`x')*(ivtheta_reg[`x',1]*(livdeltamid[1,2]*avgw`y'light+livbetamid[1,1]*ivldelta_top[1,1]*avgw`y'light_all ));
	};
};
mat ivlreg_lightelas=J(5,4,0);

forval x=1/`nbrands1'{;
	forval y=1/`nbrands2'{;
		mat ivlreg_lightelas[`x',`y']=ivlinereg`x'light`y';
	};
};



forval x=1/`nbrands2'{;
	forval y=1/`nbrands1'{;
		scalar ivlinelight`x'reg`y'=(p1`y'/w2`x')*(ivtheta_light[`x',1]*(livdeltamid[2,1]*avgw`y'reg+livbetamid[2,1]*ivldelta_top[1,1]*avgw`y'reg_all ));
	};
};

mat ivllight_regelas=J(4,5,0);

forval x=1/`nbrands2'{;
	forval y=1/`nbrands1'{;
		mat ivllight_regelas[`x',`y']=ivlinelight`x'reg`y';
	};
};

/*Linear Elasticities estimated by 2SLS*/
mat ivlinelasticities=J(9,9,0);
mat ivlinelasticities[1,1]=ivlregelas;
mat ivlinelasticities[6,6]=ivllightelas;
mat ivlinelasticities[1,6]=ivlreg_lightelas;
mat ivlinelasticities[6,1]=ivllight_regelas;


/*s_i: creates average share vector for each region*/
#delimit;
local nregions=47;
local nbrands=9;
/*Average shares by region for simulations*/
forval x=1/`nregions'{;
    mat define sbar`x'=J(`nbrands',1,0);
    forval y=1/`nbrands'{;
        qui sum share if reg`x'==1 & bra`y'==1;
        matrix sbar`x'[`y',1]=_result(3);
	};
};

/*Average expenditures by segment and region*/
mat define x_reg=J(47,1,0);
mat define x_light=J(47,1,0);

forval x=1/`nregions'{;
    qui sum expmid if reg`x'==1 & nest==1;
	matrix x_reg[`x',1]=_result(3);
	qui sum expmid if reg`x'==1 & nest==2;
	matrix x_light[`x',1]=_result(3);
};


mat X=[x_reg,x_light];


/*Revenue share weights by region*/
mat define w_nest=J(47,9,0);

forval x=1/`nregions'{;
	forval y=1/5{;
		mat w_nest[`x',`y']=w`x'b`y'reg;
	};
};
forval x=1/`nregions'{;
	forval y=1/4{;
		mat w_nest[`x',`y'+5]=w`x'b`y'light;
	};
};

mat define w=J(47,9,0);
forval x=1/`nregions'{;
	forval y=1/5{;
		mat w[`x',`y']=w`x'b`y'all;
	};
};
forval x=1/`nregions'{;
	forval y=1/4{;
		mat w[`x',`y'+5]=w`x'b`y'all;
	};
};

/*For use in bootstrapped elasticity evaluations*/
mat wnest=J(9,1,0);
forval x=1/5{;
	mat wnest[`x',1]=avgw`x'reg;
};
forval x=1/4{;
	mat wnest[5+`x',1]=avgw`x'light;
};	

mat wall=J(9,1,0);
forval x=1/5{;
	mat wall[`x',1]=avgw`x'reg_all;
};
forval x=1/4{;
	mat wall[5+`x',1]=avgw`x'light_all;
};	






/*creates average prices for each region*/
forval x=1/`nregions'{;
    mat define pbar`x'=J(`nbrands',1,0);
    forval y=1/`nbrands'{;
        qui sum price if reg`x'==1 & bra`y'==1;
        matrix pbar`x'[`y',1]=_result(3);
       };
};




/*
mat define tsale=J(`nregions',1,0);
forval x=1/`nregions'{;
    qui sum exp if reg`x'==1;
    matrix tsale[`x',1]=_result(3);
};

gen LTSALE=log(exp);
mat define ltsale=J(`nregions',1,0);
forval x=1/`nregions'{;
    qui sum LTSALE if reg`x'==1;
    matrix ltsale[`x',1]=_result(3);
};


/*creates average log prices for each region*/
forval x=1/`nregions'{;
    mat define lpbar`x'=J(`nbrands',1,0);
    forval y=1/`nbrands'{;
        qui sum lprice if bra`y'==1 & reg`x'==1;
        matrix lpbar`x'[`y',1]=_result(3);
        };
    };
*/	



/*creates average volume for each region*/
forval x=1/`nregions'{;
    mat define vbar`x'=J(`nbrands',1,0);
    forval y=1/`nbrands'{;
        qui sum wvol if bra`y'==1 & reg`x'==1;
        matrix vbar`x'[`y',1]=_result(3);
        };
    };


/*makes pre-merger "ownership" matrices*/
forval i=1/`nregions'{;
    mat define aidomega`i'=J(`nbrands',`nbrands',0);
};

forval i=1/`nregions'{;
    mat define ivaidomega`i'=J(`nbrands',`nbrands',0);
};
forval i=1/`nregions'{;
    mat define X`i'=J(`nbrands',`nbrands',0);
};

/*slope parameters for full solution*/
mat define g_reg=J(`nbrands1',`nbrands1',0);
/*[x,y] element of g_reg is the coefficient on log price of good y in equation x*/
forval x=1/5{;
	forval y=1/5{;
		mat g_reg[`x',`y']=g`x'p`y';
	};
};
mat define g_light=J(`nbrands2',`nbrands2',0);
forval x=1/4{;
	forval y=1/4{;
		mat g_light[`x',`y']=g2`x'p`y';
	};
};

mat define ivg_reg=J(`nbrands1',`nbrands1',0);
forval x=1/5{;
	forval y=1/5{;
		mat ivg_reg[`x',`y']=ivg`x'p`y';
	};
};
mat define ivg_light=J(`nbrands2',`nbrands2',0);
forval x=1/4{;
	forval y=1/4{;
		mat ivg_light[`x',`y']=ivg2`x'p`y';
	};
}; 

mat define b_reg=J(`nbrands1',1,0);
forval x=1/5{;
	mat b_reg[`x',1]=b`x';
};	

mat define b_light=J(`nbrands2',1,0);
forval x=1/4{;
	mat b_light[`x',1]=b2`x';
};	

mat define ivb_reg=J(`nbrands1',1,0);
forval x=1/5{;
	mat ivb_reg[`x',1]=ivb`x';
};	

mat define ivb_light=J(`nbrands2',1,0);
forval x=1/4{;
	mat ivb_light[`x',1]=ivb2`x';
};	
mat delta=[delta_mid_reg_reg, delta_mid_reg_light\ delta_mid_light_reg, delta_mid_light_light];
mat ivdelta=[ivdelta_mid_reg_reg, ivdelta_mid_reg_light\ ivdelta_mid_light_reg, ivdelta_mid_light_light];

mat beta=[beta_reg\beta_light];
mat ivbeta=[ivbeta_reg\ivbeta_light];

mat delta_top=[delta_top];
mat ivdelta_top=[ivdelta_top];
 
 
/*diagonal of ivaidomega matrix.*/ 
forval i=1/`nregions'{;
	forval x=1/`nbrands1'{;
			mat ivaidomega`i'[`x',`x']=(-1+(1+ivdelta_mid_reg_reg)*w_nest[`i',`x']
			+ivbeta_reg*(1+ivdelta_top)*w[`i',`x']
			+ivg`x'p`x'/sbar`i'[`x',1]
			+(ivb`x'/sbar`i'[`x',1])*(ivdelta_mid_reg_reg*w_nest[`i',`x']+w[`i',`x']*(1+ivdelta_top)*ivbeta_reg))*x_reg[`i',1]*sbar`i'[`x',1];
	};
};

forval i=1/`nregions'{;
	forval x=1/`nbrands2'{;
			mat ivaidomega`i'[5+`x',5+`x']=(-1+(1+ivdelta_mid_light_light)*w_nest[`i',5+`x']
			+ivbeta_light*(1+ivdelta_top)*w[`i',5+`x']
			+ivg2`x'p`x'/sbar`i'[5+`x',1]
			+(ivb2`x'/sbar`i'[5+`x',1])*(ivdelta_mid_light_light*w_nest[`i',5+`x']+w[`i',5+`x']*(1+ivdelta_top)*ivbeta_light))*x_light[`i',1]*sbar`i'[`x'+5,1];
	};
};



/*Need to fill out off diagonal of aidomega and ivaidomega matrices.  */
forval i=1/`nregions'{;
	mat ivaidomega`i'[6,1]=((1+ivb1/sbar`i'[1,1])*(ivbeta_reg*(1+ivdelta_top)*w[`i',6]+ivdelta_mid_reg_light*w_nest[`i',6]))*x_reg[`i',1]*sbar`i'[1,1];
	};
forval i=1/`nregions'{;
	mat ivaidomega`i'[1,6]=((1+ivb21/sbar`i'[6,1])*(ivbeta_light*(1+ivdelta_top)*w[`i',1]+ivdelta_mid_light_reg*w_nest[`i',1]))*x_light[`i',1]*sbar`i'[6,1];
};
forval i=1/`nregions'{;
	mat ivaidomega`i'[7,2]=((1+ivb2/sbar`i'[2,1])*(ivbeta_reg*(1+ivdelta_top)*w[`i',7]+ivdelta_mid_reg_light*w_nest[`i',7]))*x_reg[`i',1]*sbar`i'[2,1];
};
forval i=1/`nregions'{;
	mat ivaidomega`i'[2,7]=((1+ivb22/sbar`i'[7,1])*(ivbeta_light*(1+ivdelta_top)*w[`i',2]+ivdelta_mid_light_reg*w_nest[`i',2]))*x_light[`i',1]*sbar`i'[7,1];
};
forval i=1/`nregions'{;
	mat ivaidomega`i'[8,3]=((1+ivb3/sbar`i'[3,1])*(ivbeta_reg*(1+ivdelta_top)*w[`i',8]+ivdelta_mid_reg_light*w_nest[`i',8]))*x_reg[`i',1]*sbar`i'[3,1];
};
forval i=1/`nregions'{;
	mat ivaidomega`i'[3,8]=((1+ivb23/sbar`i'[8,1])*(ivbeta_light*(1+ivdelta_top)*w[`i',3]+ivdelta_mid_light_reg*w_nest[`i',3]))*x_light[`i',1]*sbar`i'[8,1];
};
forval i=1/`nregions'{;
	mat ivaidomega`i'[9,4]=((1+ivb4/sbar`i'[4,1])*(ivbeta_reg*(1+ivdelta_top)*w[`i',9]+ivdelta_mid_reg_light*w_nest[`i',9]))*x_reg[`i',1]*sbar`i'[4,1];
};
forval i=1/`nregions'{;
	mat ivaidomega`i'[4,9]=((1+ivb24/sbar`i'[9,1])*(ivbeta_light*(1+ivdelta_top)*w[`i',4]+ivdelta_mid_light_reg*w_nest[`i',4]))*x_light[`i',1]*sbar`i'[9,1];
};




forval i=1/`nregions'{;
	forval x=1/`nbrands1'{;
			mat aidomega`i'[`x',`x']=(-1+(1+delta_mid_reg_reg)*w_nest[`i',`x']
			+beta_reg*(1+delta_top)*w[`i',`x']
			+g`x'p`x'/sbar`i'[`x',1]
			+(b`x'/sbar`i'[`x',1])*(delta_mid_reg_reg*w_nest[`i',`x']+w[`i',`x']*(1+delta_top)*beta_reg))*x_reg[`i',1]*sbar`i'[`x',1];
	};
};

forval i=1/`nregions'{;
	forval x=1/`nbrands2'{;
			mat aidomega`i'[5+`x',5+`x']=(-1+(1+delta_mid_light_light)*w_nest[`i',5+`x']
			+beta_light*(1+delta_top)*w[`i',5+`x']
			+g2`x'p`x'/sbar`i'[5+`x',1]
			+(b2`x'/sbar`i'[5+`x',1])*(delta_mid_light_light*w_nest[`i',5+`x']+w[`i',5+`x']*(1+delta_top)*beta_light))*x_light[`i',1]*sbar`i'[`x'+5,1];
	};
};


forval i=1/`nregions'{;
	mat aidomega`i'[6,1]=((1+b1/sbar`i'[1,1])*(beta_reg*(1+delta_top)*w[`i',6]+delta_mid_reg_light*w_nest[`i',6]))*x_reg[`i',1]*sbar`i'[1,1];
	};
forval i=1/`nregions'{;
	mat aidomega`i'[1,6]=((1+b21/sbar`i'[6,1])*(beta_light*(1+delta_top)*w[`i',1]+delta_mid_light_reg*w_nest[`i',1]))*x_light[`i',1]*sbar`i'[6,1];
};
forval i=1/`nregions'{;
	mat aidomega`i'[7,2]=((1+b2/sbar`i'[2,1])*(beta_reg*(1+delta_top)*w[`i',7]+delta_mid_reg_light*w_nest[`i',7]))*x_reg[`i',1]*sbar`i'[2,1];
};
forval i=1/`nregions'{;
	mat aidomega`i'[2,7]=((1+b22/sbar`i'[7,1])*(beta_light*(1+delta_top)*w[`i',2]+delta_mid_light_reg*w_nest[`i',2]))*x_light[`i',1]*sbar`i'[7,1];
};
forval i=1/`nregions'{;
	mat aidomega`i'[8,3]=((1+b3/sbar`i'[3,1])*(beta_reg*(1+delta_top)*w[`i',8]+delta_mid_reg_light*w_nest[`i',8]))*x_reg[`i',1]*sbar`i'[3,1];
};
forval i=1/`nregions'{;
	mat aidomega`i'[3,8]=((1+b23/sbar`i'[8,1])*(beta_light*(1+delta_top)*w[`i',3]+delta_mid_light_reg*w_nest[`i',3]))*x_light[`i',1]*sbar`i'[8,1];
};
forval i=1/`nregions'{;
	mat aidomega`i'[9,4]=((1+b4/sbar`i'[4,1])*(beta_reg*(1+delta_top)*w[`i',9]+delta_mid_reg_light*w_nest[`i',9]))*x_reg[`i',1]*sbar`i'[4,1];
};
forval i=1/`nregions'{;
	mat aidomega`i'[4,9]=((1+b24/sbar`i'[9,1])*(beta_light*(1+delta_top)*w[`i',4]+delta_mid_light_reg*w_nest[`i',4]))*x_light[`i',1]*sbar`i'[9,1];
};

forval i=1/`nregions'{;
    forval x=1/`nbrands1'{;
        mat X`i'[`x',`x']=x_reg[`i',1];
    };
};

forval i=1/`nregions'{;
    forval x=1/`nbrands2'{;
        mat X`i'[`x'+5,`x'+5]=x_light[`i',1];
    };
};




/*aidmu has p-mc, not lerner index*/
forval i=1/`nregions'{;
    mat aidmu`i'=J(`nbrands',1,0);
    mat aidmu`i'=-inv(aidomega`i')*X`i'*sbar`i';
    mat ivaidmu`i'=J(`nbrands',1,0);
    mat ivaidmu`i'=-inv(ivaidomega`i')*X`i'*sbar`i';
};



/*mark-up matrix has 10 columns, which correspond to regional markets, and 7 rows corresponding to brands.  
aidmu are p-mc's with aid system, lmu are p=mc's with linear system, llmu are lerner indices with log-linear system
this just puts them in matrix form as described above*/
mat aidmu=[aidmu1];
mat ivaidmu=[ivaidmu1];
mat pbar=[pbar1];
/*mat lpbar=[lpbar1];*/
mat sbar=[sbar1];
forval x=2/`nregions'{;
    mat aidmu=[aidmu,aidmu`x'];
    mat ivaidmu=[ivaidmu,ivaidmu`x'];
    mat pbar=[pbar,pbar`x'];
    mat sbar=[sbar,sbar`x'];
};

mat pbar=pbar';
mat sbar=sbar';
mat aidmu=aidmu';
mat ivaidmu=ivaidmu';

/*get costs by solving mark-ups.  because element-by-element operations are a pain in stata, i do this in mata*/
/*input markup matrices, average prices and shares, and post-merger ownership matrices into stata*/
mata;
	b_reg=st_matrix("b_reg");
	b_light=st_matrix("b_light");
	ivb_reg=st_matrix("ivb_reg");
	ivb_light=st_matrix("ivb_light");
	
	g_reg=st_matrix("g_reg");
	g_light=st_matrix("g_light");
	ivg_reg=st_matrix("ivg_reg");
	ivg_light=st_matrix("ivg_light");
	
	w=st_matrix("w");
	w_nest=st_matrix("w_nest");
	
	X=st_matrix("X");
	
	delta=st_matrix("delta");
	ivdelta=st_matrix("ivdelta");
	
	beta=st_matrix("beta");
	ivbeta=st_matrix("ivbeta");
	
	delta_top=st_matrix("delta_top");
	ivdelta_top=st_matrix("ivdelta_top");
	
    aidmu=st_matrix("aidmu");   
    ivaidmu=st_matrix("ivaidmu");

    pbar=st_matrix("pbar");
    sbar=st_matrix("sbar");
    aidcost=pbar:*(J(`nregions',`nbrands',1)-aidmu);
    ivaidcost=pbar:*(J(`nregions',`nbrands',1)-ivaidmu);
    st_matrix("aidcost", aidcost);
    st_matrix("ivaidcost", ivaidcost);
end;



/*defines post-merger ownership matrix of elasticities*/
forval i=1/`nregions'{;
    mat define taidomega`i'=J(`nbrands',`nbrands',0);
};
forval i=1/`nregions'{;
    mat define ivtaidomega`i'=J(`nbrands',`nbrands',0);
};


forval i=1/`nregions'{;
	mat ivtaidomega`i'[8,4]=((1+ivb4/sbar`i'[4,1])*(ivbeta_reg*(1+ivdelta_top)*w[`i',8]+ivdelta_mid_reg_light*w_nest[`i',8]))*x_reg[`i',1]*sbar`i'[4,1];
	};
forval i=1/`nregions'{;
	mat ivtaidomega`i'[4,8]=((1+ivb23/sbar`i'[8,1])*(ivbeta_light*(1+ivdelta_top)*w[`i',4]+ivdelta_mid_light_reg*w_nest[`i',4]))*x_light[`i',1]*sbar`i'[8,1];
};

forval i=1/`nregions'{;
	mat ivtaidomega`i'[9,3]=((1+ivb3/sbar`i'[3,1])*(ivbeta_reg*(1+ivdelta_top)*w[`i',9]+ivdelta_mid_reg_light*w_nest[`i',9]))*x_reg[`i',1]*sbar`i'[3,1];
};
forval i=1/`nregions'{;
	mat ivtaidomega`i'[3,9]=((1+ivb24/sbar`i'[9,1])*(ivbeta_light*(1+ivdelta_top)*w[`i',3]+ivdelta_mid_light_reg*w_nest[`i',3]))*x_light[`i',1]*sbar`i'[9,1];
};



forval i=1/`nregions'{;            
			mat ivtaidomega`i'[4,3]=((1+ivdelta_mid_reg_reg)*w_nest[`i',4]
			+ivbeta_reg*(1+ivdelta_top)*w[`i',4]
			+ivg3p4/sbar`i'[3,1]
			+(ivb3/sbar`i'[3,1])*(ivdelta_mid_reg_reg*w_nest[`i',4]+w[`i',4]*(1+ivdelta_top)*ivbeta_reg))*x_reg[`i',1]*sbar`i'[3,1];
	};
	

forval i=1/`nregions'{;
			mat ivtaidomega`i'[3,4]=((1+ivdelta_mid_reg_reg)*w_nest[`i',3]
			+ivbeta_reg*(1+ivdelta_top)*w[`i',3]
			+ivg4p3/sbar`i'[4,1]
			+(ivb4/sbar`i'[4,1])*(ivdelta_mid_reg_reg*w_nest[`i',3]+w[`i',3]*(1+ivdelta_top)*ivbeta_reg))*x_reg[`i',1]*sbar`i'[4,1];
};


forval i=1/`nregions'{;
			mat ivtaidomega`i'[9,8]=((1+ivdelta_mid_light_light)*w_nest[`i',9]
			+ivbeta_light*(1+ivdelta_top)*w[`i',9]
			+ivg23p4/sbar`i'[8,1]
			+(ivb23/sbar`i'[8,1])*(ivdelta_mid_light_light*w_nest[`i',9]+w[`i',9]*(1+ivdelta_top)*ivbeta_light))*x_light[`i',1]*sbar`i'[8,1];
};

forval i=1/`nregions'{;
	mat ivtaidomega`i'[8,9]=((1+ivdelta_mid_light_light)*w_nest[`i',8]
			+ivbeta_light*(1+ivdelta_top)*w[`i',8]
			+ivg24p3/sbar`i'[9,1]
			+(ivb24/sbar`i'[9,1])*(ivdelta_mid_light_light*w_nest[`i',8]+w[`i',8]*(1+ivdelta_top)*ivbeta_light))*x_light[`i',1]*sbar`i'[9,1];
};


forval i=1/`nregions'{;
	mat taidomega`i'[8,4]=((1+b4/sbar`i'[4,1])*(beta_reg*(1+delta_top)*w[`i',8]+delta_mid_reg_light*w_nest[`i',8]))*x_reg[`i',1]*sbar`i'[4,1];
	};


forval i=1/`nregions'{;
	mat taidomega`i'[4,8]=((1+b23/sbar`i'[8,1])*(beta_light*(1+delta_top)*w[`i',4]+delta_mid_light_reg*w_nest[`i',4]))*x_light[`i',1]*sbar`i'[8,1];
};

forval i=1/`nregions'{;
	mat taidomega`i'[9,3]=((1+b3/sbar`i'[3,1])*(beta_reg*(1+delta_top)*w[`i',9]+delta_mid_reg_light*w_nest[`i',9]))*x_reg[`i',1]*sbar`i'[3,1];
};


forval i=1/`nregions'{;
	mat taidomega`i'[3,9]=((1+b24/sbar`i'[9,1])*(beta_light*(1+delta_top)*w[`i',3]+delta_mid_light_reg*w_nest[`i',3]))*x_light[`i',1]*sbar`i'[9,1];
};



forval i=1/`nregions'{;            
			mat taidomega`i'[4,3]=((1+delta_mid_reg_reg)*w_nest[`i',4]
			+beta_reg*(1+delta_top)*w[`i',4]
			+g3p4/sbar`i'[3,1]
			+(b3/sbar`i'[3,1])*(delta_mid_reg_reg*w_nest[`i',4]+w[`i',4]*(1+delta_top)*beta_reg))*x_reg[`i',1]*sbar`i'[3,1];
	};

	
forval i=1/`nregions'{;
			mat taidomega`i'[3,4]=((1+delta_mid_reg_reg)*w_nest[`i',3]
			+beta_reg*(1+delta_top)*w[`i',3]
			+g4p3/sbar`i'[4,1]
			+(b4/sbar`i'[4,1])*(delta_mid_reg_reg*w_nest[`i',3]+w[`i',3]*(1+delta_top)*beta_reg))*x_reg[`i',1]*sbar`i'[4,1];
};



forval i=1/`nregions'{;
			mat taidomega`i'[9,8]=((1+delta_mid_light_light)*w_nest[`i',9]
			+beta_light*(1+delta_top)*w[`i',9]
			+g23p4/sbar`i'[8,1]
			+(b23/sbar`i'[8,1])*(delta_mid_light_light*w_nest[`i',9]+w[`i',9]*(1+delta_top)*beta_light))*x_light[`i',1]*sbar`i'[8,1];
};



forval i=1/`nregions'{;
	mat taidomega`i'[8,9]=((1+delta_mid_light_light)*w_nest[`i',8]
			+beta_light*(1+delta_top)*w[`i',8]
			+g24p3/sbar`i'[9,1]
			+(b24/sbar`i'[9,1])*(delta_mid_light_light*w_nest[`i',8]+w[`i',8]*(1+delta_top)*beta_light))*x_light[`i',1]*sbar`i'[9,1];
};





#delimit;
/*solves for post-merger mark-ups*/
forval i=1/`nregions'{;
    mat postaidomega`i'=aidomega`i'+taidomega`i';
};
forval i=1/`nregions'{;
    mat ivpostaidomega`i'=ivaidomega`i'+ivtaidomega`i';
};

forval i=1/`nregions'{;
    mat paidmu`i'=-inv(postaidomega`i')*X`i'*sbar`i';
};

forval i=1/`nregions'{;
    mat ivpaidmu`i'=-inv(ivpostaidomega`i')*X`i'*sbar`i';
};

/*puts post-merger mark-ups into a matrix with rows corresponding to markets and columns to products*/

mat paidmu=[paidmu1];
mat ivpaidmu=[ivpaidmu1];

forval x=2/`nregions'{;
    mat paidmu=[paidmu,paidmu`x'];
};

mat paidmu=paidmu';


forval x=2/`nregions'{;
    mat ivpaidmu=[ivpaidmu,ivpaidmu`x'];
};

mat ivpaidmu=ivpaidmu';


# delimit;



/*now go back to mata to solve for post-merger prices.  will then calculate percentage increases over 
pre-merger average prices for each region, and then take the average over regions for to yield percentage
price increases for each product.  these average prices are the mxprice variables and the average costs are the
mxcost variables*/

mata;
    /*tsale=st_matrix("tsale");*/
    paidmu=st_matrix("paidmu");
    ivpaidmu=st_matrix("ivpaidmu");
    invaidp=(J(`nregions',`nbrands',1)-paidmu):/aidcost;
    ivinvaidp=(J(`nregions',`nbrands',1)-ivpaidmu):/ivaidcost;
    aidp=((J(`nregions',`nbrands',1):/invaidp)-pbar):/pbar;
    ivaidp=((J(`nregions',`nbrands',1):/ivinvaidp)-pbar):/pbar;
    meanaidcost=mean(aidcost,1);
    ivmeanaidcost=mean(ivaidcost,1);
    meanaidprice=mean(aidp,1);
    ivmeanaidprice=mean(ivaidp,1);
    medianaidprice=mm_median(aidp);
    ivmedianaidprice=mm_median(ivaidp);
	postdirect=J(47,9,0)
	postdirect[.,1]=pbar[.,1]*(1+0.0003)  /*aj*/
	postdirect[.,2]=pbar[.,2]*(1+0.0160) /*hj*/
	postdirect[.,3]=pbar[.,3]*(1+0.0323)  /*lj*/
	postdirect[.,4]=pbar[.,4]*(1-0.0090)  /*mv*/
	postdirect[.,5]=pbar[.,5]*(1+0.0072)  /*pl*/
	postdirect[.,6]=pbar[.,6]*(1+0.0094)  /*ajl*/
	postdirect[.,7]=pbar[.,7]*(1+0.0172)  /*hjl*/
	postdirect[.,8]=pbar[.,8]*(1+0.0346)  /*lvl*/
	postdirect[.,9]=pbar[.,9]*(1+0.0037)  /*mbl*/
end;    



/****************SOLVES UNDER LINEAR DEMAND*****************************/


#delimit cr
/*mata matuse c:/data/restat_12_9/analysis_data/post_prices, replace*/
mata:
		
    B_reg=st_matrix("B_reg")
	B_light=st_matrix("B_light")
    V=(st_matrix("vbar1"),st_matrix("vbar2"),st_matrix("vbar3"),st_matrix("vbar4"),st_matrix("vbar5"),st_matrix("vbar6"),st_matrix("vbar7"),st_matrix("vbar8"),st_matrix("vbar9"),st_matrix("vbar10"),
	st_matrix("vbar11"),st_matrix("vbar12"),st_matrix("vbar13"),st_matrix("vbar14"),st_matrix("vbar15"),st_matrix("vbar16"),st_matrix("vbar17"),st_matrix("vbar18"),st_matrix("vbar19"),st_matrix("vbar20"),
	st_matrix("vbar21"),st_matrix("vbar22"),st_matrix("vbar23"),st_matrix("vbar24"),st_matrix("vbar25"),st_matrix("vbar26"),st_matrix("vbar27"),st_matrix("vbar28"),st_matrix("vbar29"),st_matrix("vbar30"),
	st_matrix("vbar31"),st_matrix("vbar32"),st_matrix("vbar33"),st_matrix("vbar34"),st_matrix("vbar35"),st_matrix("vbar36"),st_matrix("vbar37"),st_matrix("vbar38"),st_matrix("vbar39"),st_matrix("vbar40"),
	st_matrix("vbar41"),st_matrix("vbar42"),st_matrix("vbar43"),st_matrix("vbar44"),st_matrix("vbar45"),st_matrix("vbar46"),st_matrix("vbar47"))
	theta_reg=st_matrix("theta_reg")
	theta_light=st_matrix("theta_light")
	ldelta_mid=st_matrix("ldeltamid")
	lbeta_mid=st_matrix("lbetamid")
	ldelta_top=st_matrix("ldelta_top")
	lbeta_reg=st_matrix("lbeta_reg")
		
	
	a=J(9,47,0) 
    c=J(9,47,0)
    pbar=pbar'
	postdirect=postdirect'
	/*create intercepts for regular and then light*/
    for(i=1;i<=47;i++){
        a[1..5,i]=V[1..5,i]-B_reg*pbar[1..5,i]       
    }
	for(i=1;i<=47;i++){
        a[6..9,i]=V[6..9,i]-B_light*pbar[6..9,i]       
    }
	
	B=J(9,9,0)
	B[1..5,1..5]=B_reg
	B[6..9,6..9]=B_light
	/*Create slope matrix here, accounting for middle and upper levels of demand*/
	
	
	mc=J(47,9,0)
	mcpost=J(47,9,0)
	lfullp=J(47,9,0)
	for(j=1;j<=47;j++){
		D_pre=J(9,9,0)
		D_post=J(9,9,0)
		temp=J(9,9,0)
		for(i=1;i<=5;i++){
			D_pre[i,i]=B_reg[i,i]+theta_reg[i,1]*(ldelta_mid[1,1]*w_nest[j,i]+w[j,i]*ldelta_top*lbeta_mid[1,1])
		}
		for(i=6;i<=9;i++){
			D_pre[i,i]=B_light[i-5,i-5]+theta_light[i-5,1]*(ldelta_mid[2,2]*w_nest[j,i]+w[j,i]*ldelta_top*lbeta_mid[2,1])
		}
		D_pre[1,6]=theta_light[1,1]*(ldelta_mid[2,1]*w_nest[j,1]+w[j,1]*ldelta_top*lbeta_mid[2,1])
		D_pre[6,1]=theta_reg[1,1]*(ldelta_mid[1,2]*w_nest[j,6]+w[j,6]*ldelta_top*lbeta_mid[1,1])
		D_pre[2,7]=theta_light[2,1]*(ldelta_mid[2,1]*w_nest[j,2]+w[j,2]*ldelta_top*lbeta_mid[2,1])
		D_pre[7,2]=theta_reg[2,1]*(ldelta_mid[1,2]*w_nest[j,7]+w[j,7]*ldelta_top*lbeta_mid[1,1])
		D_pre[3,8]=theta_light[3,1]*(ldelta_mid[2,1]*w_nest[j,3]+w[j,3]*ldelta_top*lbeta_mid[2,1])
		D_pre[8,3]=theta_reg[3,1]*(ldelta_mid[1,2]*w_nest[j,8]+w[j,8]*ldelta_top*lbeta_mid[1,1])
		D_pre[4,9]=theta_light[4,1]*(ldelta_mid[2,1]*w_nest[j,4]+w[j,4]*ldelta_top*lbeta_mid[2,1])
		D_pre[9,4]=theta_reg[4,1]*(ldelta_mid[1,2]*w_nest[j,9]+w[j,9]*ldelta_top*lbeta_mid[1,1])
		
		temp[3,4]=B_reg[4,3]+theta_reg[4,1]*(ldelta_mid[1,1]*w_nest[j,3]+w[j,3]*ldelta_top*lbeta_mid[1,1])
		temp[4,3]=B_reg[3,4]+theta_reg[3,1]*(ldelta_mid[1,1]*w_nest[j,4]+w[j,4]*ldelta_top*lbeta_mid[1,1])
		temp[8,9]=B_light[4,3]+theta_light[4,1]*(ldelta_mid[2,2]*w_nest[j,8]+w[j,8]*ldelta_top*lbeta_mid[2,1])
		temp[9,8]=B_light[3,4]+theta_light[3,1]*(ldelta_mid[2,2]*w_nest[j,9]+w[j,9]*ldelta_top*lbeta_mid[2,1])
		
		temp[3,9]=theta_light[3,1]*(ldelta_mid[2,1]*w_nest[j,3]+w[j,3]*ldelta_top*lbeta_mid[2,1])
		temp[9,3]=theta_reg[3,1]*(ldelta_mid[1,2]*w_nest[j,9]+w[j,9]*ldelta_top*lbeta_mid[1,1])
		temp[4,8]=theta_light[4,1]*(ldelta_mid[2,1]*w_nest[j,4]+w[j,4]*ldelta_top*lbeta_mid[2,1])
		temp[8,4]=theta_reg[4,1]*(ldelta_mid[1,2]*w_nest[j,8]+w[j,8]*ldelta_top*lbeta_mid[1,1])
		
		D_post=temp+D_pre
		
		mcind=luinv(D_pre)*(a[.,j]+B*pbar[.,j]+D_pre*pbar[.,j])
		mc[j,.]=mcind'
		pind=luinv(B+D_post)*(D_post*mc[j,.]'-a[.,j])
		lfullp[j,.]=pind'		
		mcpostind=luinv(D_post)*(a[.,j]+B*postdirect[.,j]+D_post*postdirect[.,j])
		mcpost[j,.]=mcpostind'
	}	
	linans=mm_median((lfullp-pbar'):/pbar')
	mcchange=mm_median((mcpost-mc):/mc)
end



mata:
	ivmcpost=J(47,9,0)
    ivB_reg=st_matrix("ivB_reg")
	ivB_light=st_matrix("ivB_light")    
	ivtheta_reg=st_matrix("ivtheta_reg")
	ivtheta_light=st_matrix("ivtheta_light")
	ivldelta_mid=st_matrix("livdeltamid")
	ivlbeta_mid=st_matrix("livbetamid")
	ivldelta_top=st_matrix("ivldelta_top")
	ivlbeta_reg=st_matrix("livbeta_reg")
	a=J(9,47,0) 
    c=J(9,47,0)
    
	/*create intercepts for regular and then light*/
    for(i=1;i<=47;i++){
        a[1..5,i]=V[1..5,i]-ivB_reg*pbar[1..5,i]       
    }
	for(i=1;i<=47;i++){
        a[6..9,i]=V[6..9,i]-ivB_light*pbar[6..9,i]       
    }
	ivB=J(9,9,0)
	ivB[1..5,1..5]=ivB_reg
	ivB[6..9,6..9]=ivB_light
	/*Create slope matrix here, accounting for middle and upper levels of demand*/
	ivmc=J(47,9,0)
	ivlfullp=J(47,9,0)
	for(j=1;j<=47;j++){
		ivD_pre=J(9,9,0)
		ivD_post=J(9,9,0)
		temp=J(9,9,0)
		for(i=1;i<=5;i++){
			ivD_pre[i,i]=ivB_reg[i,i]+ivtheta_reg[i,1]*(ivldelta_mid[1,1]*w_nest[j,i]+w[j,i]*ivldelta_top*ivlbeta_mid[1,1])
		}
		for(i=6;i<=9;i++){
			ivD_pre[i,i]=ivB_light[i-5,i-5]+ivtheta_light[i-5,1]*(ivldelta_mid[2,2]*w_nest[j,i]+w[j,i]*ivldelta_top*ivlbeta_mid[2,1])
		}	
		
		ivD_pre[1,6]=ivtheta_light[1,1]*(ivldelta_mid[2,1]*w_nest[j,1]+w[j,1]*ivldelta_top*ivlbeta_mid[2,1])
		ivD_pre[6,1]=ivtheta_reg[1,1]*(ivldelta_mid[1, 2]*w_nest[j,6]+w[j,6]*ivldelta_top*ivlbeta_mid[1,1])
		ivD_pre[2,7]=ivtheta_light[2,1]*(ivldelta_mid[2,1]*w_nest[j,2]+w[j,2]*ivldelta_top*ivlbeta_mid[2,1])
		ivD_pre[7,2]=ivtheta_reg[2,1]*(ivldelta_mid[1,2]*w_nest[j,7]+w[j,7]*ivldelta_top*ivlbeta_mid[1,1])
		ivD_pre[3,8]=ivtheta_light[3,1]*(ivldelta_mid[2,1]*w_nest[j,3]+w[j,3]*ivldelta_top*ivlbeta_mid[2,1])
		ivD_pre[8,3]=ivtheta_reg[3,1]*(ivldelta_mid[1,2]*w_nest[j,8]+w[j,8]*ivldelta_top*ivlbeta_mid[1,1])
		ivD_pre[4,9]=ivtheta_light[4,1]*(ivldelta_mid[2,1]*w_nest[j,4]+w[j,4]*ivldelta_top*ivlbeta_mid[2,1])
		ivD_pre[9,4]=ivtheta_reg[4,1]*(ivldelta_mid[1,2]*w_nest[j,9]+w[j,9]*ivldelta_top*ivlbeta_mid[1,1])
		
		temp[3,4]=ivB_reg[4,3]+ivtheta_reg[4,1]*(ivldelta_mid[1,1]*w_nest[j,3]+w[j,3]*ivldelta_top*ivlbeta_mid[1,1])
		temp[4,3]=ivB_reg[3,4]+ivtheta_reg[3,1]*(ivldelta_mid[1,1]*w_nest[j,4]+w[j,4]*ivldelta_top*ivlbeta_mid[1,1])
		temp[8,9]=ivB_light[4,3]+ivtheta_light[4,1]*(ivldelta_mid[2,2]*w_nest[j,8]+w[j,8]*ivldelta_top*ivlbeta_mid[2,1])
		temp[9,8]=ivB_light[3,4]+ivtheta_light[3,1]*(ivldelta_mid[2,2]*w_nest[j,9]+w[j,9]*ivldelta_top*ivlbeta_mid[2,1])
		
		temp[3,9]=ivtheta_light[3,1]*(ivldelta_mid[2,1]*w_nest[j,3]+w[j,3]*ivldelta_top*ivlbeta_mid[2,1])
		temp[9,3]=ivtheta_reg[3,1]*(ivldelta_mid[1,2]*w_nest[j,9]+w[j,9]*ivldelta_top*ivlbeta_mid[1,1])
		temp[4,8]=ivtheta_light[4,1]*(ivldelta_mid[2,1]*w_nest[j,4]+w[j,4]*ivldelta_top*ivlbeta_mid[2,1])
		temp[8,4]=ivtheta_reg[4,1]*(ivldelta_mid[1,2]*w_nest[j,8]+w[j,8]*ivldelta_top*ivlbeta_mid[1,1])
		
		ivD_post=temp+ivD_pre
			
		ivmcind=luinv(ivD_pre)*(a[.,j]+ivB*pbar[.,j]+ivD_pre*pbar[.,j])
		ivmc[j,.]=ivmcind'
		ivpind=luinv(ivB+ivD_post)*(ivD_post*ivmc[j,.]'-a[.,j])
		ivlfullp[j,.]=ivpind'
		mcpostind=luinv(ivD_post)*(a[.,j]+ivB*postdirect[.,j]+ivD_post*postdirect[.,j])
		ivmcpost[j,.]=mcpostind'		
	}	
	ivlinans=mm_median((ivlfullp-pbar'):/pbar') 
	ivmcchange=mm_median((ivmcpost-ivmc):/ivmc)
end




/*this solves for the exact equilibrium for each model*/
#delimit cr
mata:
	aidfocs=J(`nregions',`nbrands',0)
	pbar=pbar'
	aidfullp=J(`nregions',`nbrands',0)
	check=J(`nregions',1,0)
    for(i=1;i<=`nregions';i++){
		a=csolve(&msyraids(),pbar[i,.]',1e-6,400,sbar[i,.]',pbar[i,.]',aidcost[i,.]',g_reg,g_light,b_reg,b_light,w[i,.]',w_nest[i,.]',X[i,.]',delta,beta,delta_top)		
        aidfullp[i,.]=a[1..`nbrands',.]'
		ws=msyraids(a[1..9,.],sbar[i,.]',pbar[i,.]',aidcost[i,.]',g_reg,g_light,b_reg,b_light,w[i,.]',w_nest[i,.]',X[i,.]',delta,beta,delta_top)
		check[i,1]=((a[10,1]+(a[1..rows(a)-1,.]==J(rows(a)-1,1,.))+(ws[1,.]==.)+(ws[2,.]==.)+(ws[3,.]==.)+(ws[4,.]==.)+(ws[5,.]==.)+(ws[6,.]==.)+(ws[7,.]==.)+(ws[8,.]==.)+(ws[9,.]==.)  )==0)
		/*aidfocs[i,.]=msyraids(postpbar[i,.]',sbar[i,.]',pbar[i,.]',aidcost[i,.]',g_reg,g_light,b_reg,b_light,w[i,.]',w_nest[i,.]',X[i,.]',delta,beta,delta_top)'*/
    }
    ans=(aidfullp-pbar):/pbar
    faids=mm_median(select(ans,check))    	
end




mata:
	ivaidfocs=J(`nregions',`nbrands',0)
	ivaidfullp=J(`nregions',`nbrands',0)
	ivcheck=J(`nregions',1,0)
    for(i=1;i<=`nregions';i++){
		a=csolve(&msyraids(),pbar[i,.]',1e-6,400,sbar[i,.]',pbar[i,.]',ivaidcost[i,.]',ivg_reg,ivg_light,ivb_reg,ivb_light,w[i,.]',w_nest[i,.]',X[i,.]',ivdelta,ivbeta,ivdelta_top)		
        ivaidfullp[i,.]=a[1..`nbrands',.]'
		ws=msyraids(a[1..9,.],sbar[i,.]',pbar[i,.]',ivaidcost[i,.]',ivg_reg,ivg_light,ivb_reg,ivb_light,w[i,.]',w_nest[i,.]',X[i,.]',ivdelta,ivbeta,ivdelta_top)
		ivcheck[i,1]=((a[10,1]+(a[1..rows(a)-1,.]==J(rows(a)-1,1,.))+(ws[1,.]==.)+(ws[2,.]==.)+(ws[3,.]==.)+(ws[4,.]==.)+(ws[5,.]==.)+(ws[6,.]==.)+(ws[7,.]==.)+(ws[8,.]==.)+(ws[9,.]==.)  )==0)
		/*ivaidfocs[i,.]=msyraids(postpbar[i,.]',sbar[i,.]',pbar[i,.]',ivaidcost[i,.]',ivg_reg,ivg_light,ivb_reg,ivb_light,w[i,.]',w_nest[i,.]',X[i,.]',ivdelta,ivbeta,ivdelta_top)'*/
    }
    ivans=(ivaidfullp-pbar):/pbar
    ivfaids=mm_median(select(ivans,ivcheck))    	
end

mata:
	linans
	ivlinans
	faids
	ivfaids
end:	





do c:/data/restat_12_9/programs/table2/mcosts_syraids
local nbrands=9
local nregions=47
mata:
    caidfullp=J(`nbrands',`nregions',0)
    for(i=1;i<=cols(pbar');i++){
		init=aidcost[i,.]'
		ca=csolve(&cost_msyraids(),init,1e-6,400,sbar[i,.]',pbar[i,.]',postdirect[.,i],g_reg,g_light,b_reg,b_light,w[i,.]',w_nest[i,.]',X[i,.]',delta,beta,delta_top)		       	
        ca=ca[1..`nbrands',.]
        caidfullp[.,i]=ca
    }
	caidfullp=caidfullp'
	mm_median(select((caidfullp[.,3]-aidcost[.,3]):/aidcost[.,3],check))
	mm_median(select((caidfullp[.,4]-aidcost[.,4]):/aidcost[.,4],check))
	mm_median(select((caidfullp[.,8]-aidcost[.,8]):/aidcost[.,8],check))
	mm_median(select((caidfullp[.,9]-aidcost[.,9]):/aidcost[.,9],check))
	mean((caidfullp[.,3]-aidcost[.,3]):/aidcost[.,3],check)     	
	mean((caidfullp[.,4]-aidcost[.,4]):/aidcost[.,4],check)    	
	mean((caidfullp[.,8]-aidcost[.,8]):/aidcost[.,8],check)
	mean((caidfullp[.,9]-aidcost[.,9]):/aidcost[.,9],check)    	
end




mata:
    ivcaidfullp=J(`nbrands',`nregions',0)
    for(i=1;i<=cols(pbar');i++){        
		init=ivaidcost[i,.]'
		ivca=csolve(&cost_msyraids(),init,1e-6,400,sbar[i,.]',pbar[i,.]',postdirect[.,i],ivg_reg,ivg_light,ivb_reg,ivb_light,w[i,.]',w_nest[i,.]',X[i,.]',ivdelta,ivbeta,ivdelta_top)		       	
        ivca=ivca[1..`nbrands',.]
        ivcaidfullp[.,i]=ivca
    }
	ivcaidfullp=ivcaidfullp'
	mm_median(select((ivcaidfullp[.,3]-ivaidcost[.,3]):/ivaidcost[.,3],ivcheck))
	mm_median(select((ivcaidfullp[.,4]-ivaidcost[.,4]):/ivaidcost[.,4],ivcheck))
	mm_median(select((ivcaidfullp[.,8]-ivaidcost[.,8]):/ivaidcost[.,8],ivcheck))
	mm_median(select((ivcaidfullp[.,9]-ivaidcost[.,9]):/ivaidcost[.,9],ivcheck))
	
	mean((ivcaidfullp[.,3]-ivaidcost[.,3]):/ivaidcost[.,3],ivcheck)
	mean((ivcaidfullp[.,4]-ivaidcost[.,4]):/ivaidcost[.,4],ivcheck)    	
	mean((ivcaidfullp[.,8]-ivaidcost[.,8]):/ivaidcost[.,8],ivcheck)
	mean((ivcaidfullp[.,9]-ivaidcost[.,9]):/ivaidcost[.,9],ivcheck)    	
end






#delimit;
save syrstart, replace; 
drop _all;

#delimit;
/***Now that we have our point estimates of costs and approximate simulated price changes, we draw from asymptotic approx of distributions
of estimators and recalculate for each draw.  will look at empirical distribution of these calculated costs and changes to do inference.****/
drop _all;
/*number of bootstrap iterations*/
local its=1000;
local aids_reg "laidpreg1 laidpreg2 laidpreg3 laidpreg4 `reglist1' `Mlist1' `regionxbrand1' bra11 bra12 bra13 bra14 count11 count12 count13 count14";
local aids_light "laidplight1 laidplight2 laidplight3 `reglist2' `Mlist2' `regionxbrand2' bra21 bra22 bra23 count21 count22 count23";
local aids_mid "`femid' `Mlistmid' nest1 nest2 lexp1 lexp2 lp1_1 lp1_2 lp2_1 lp2_2 count1 count2";
local aids_top "indexall month2 month3 month4 month5 month6 month7 month10 month11 month12 count reg2 reg3 reg4 reg5 reg6 reg7 reg8 reg9 reg10 reg11 reg12 reg13 reg14 reg15 reg16 reg17 reg18 reg19 reg20 reg21 reg22 reg23 reg24 reg25 reg26 reg27 reg28 reg29 reg30 reg31 reg32 reg33 reg34 reg35 reg36 reg37 reg38 reg39 reg40 reg41 reg42 reg43 reg44 reg45 reg46 reg47";

 

drawnorm `aids_reg', n(`its') mean(aidB) cov(aidVC);
drawnorm `aids_light', n(`its') mean(aidB2) cov(aidVC2); 
drawnorm `aids_mid', n(`its') mean(aidMidB) cov(aidMidVC);
drawnorm `aids_top', n(`its') mean(ols_b_top) cov(ols_vc_top);


#delimit;
foreach x of local list1{;
    gen lpreg_5_`x'=-lpreg_1_`x'-lpreg_2_`x'-lpreg_3_`x'-lpreg_4_`x';
};
gen laidpreg5=-laidpreg1-laidpreg2-laidpreg3-laidpreg4;

forval x=1/5{;
	rename lpreg_`x'_AUNT_JEMIMA lpreg_`x'_1;
	rename lpreg_`x'_HUNGRY_JACK lpreg_`x'_2;
	rename lpreg_`x'_LOG_CABIN lpreg_`x'_3;
	rename lpreg_`x'_MRS_BUTTERWORTH lpreg_`x'_4;
	rename lpreg_`x'_PRIVATE_LABEL lpreg_`x'_5;
};

foreach x of local list2{;
    gen lplight_4_`x'=-lplight_1_`x'-lplight_2_`x'-lplight_3_`x';
};
gen laidplight4=-laidplight1-laidplight2-laidplight3;


forval x=1/4{;
	rename lplight_`x'_AUNT_JEMIMA_LIGHT lplight_`x'_1;
	rename lplight_`x'_JACK_LITE lplight_`x'_2;
	rename lplight_`x'_LOG_CABIN_LIGHT lplight_`x'_3;
	rename lplight_`x'_MRS_B_LITE lplight_`x'_4;
};

drop M11b1-count14 M21b1-nest2 count1 count2  month2-reg47;

mkmat _all;

use syrstart, replace;
keep if merged==0;


/*bsbarx is a 7x`its' matrix with columns equal to the means of revenue shares taken over bootstrap subsamples
of region x.  bpbarx, bvbbarx defined similarly.*/

forval x=1/`nregions'{;
    mat bsbar`x'=J(9,`its',0);
    mat bpbar`x'=J(9,`its',0);
	mat x_reg`x'=J(1,`its',0);
	mat x_light`x'=J(1,`its',0);
};

forval x=1/`nregions'{;
    forval w=1/`its'{;
        preserve;
            keep if reg`x'==1;
            bsample;
			qui sum expmid if reg`x'==1&nest==1;
			mat x_reg`x'[1,`w']=_result(3);
			qui sum expmid if reg`x'==1&nest==2;
			mat x_light`x'[1,`w']=_result(3);
			
			forval y=1/`nbrands'{;
                qui sum share if bra`y'==1;
                matrix bsbar`x'[`y',`w']=_result(3);
                qui sum price if bra`y'==1;
                matrix bpbar`x'[`y',`w']=_result(3);
            };
        restore;
    };
};

#delimit;
local its=1000;
local nbrands=9;
local nregions=47;
local nbrands1=5;
local nbrands2=4;
mata;
    bsaidcost=J(`its',9,0);
    bsaidp=J(`its',9,0);
end;

#delimit;
local its=1000;
forval z=1/`its'{;
    /*makes pre-merger "ownership" matrices*/
	forval i=1/`nregions'{;
        mat define aidomega`i'=J(`nbrands',`nbrands',0);
    };

	
	
	forval i=1/`nregions'{;
		mat define X`i'=J(`nbrands',`nbrands',0);
	};

	
	forval i=1/`nregions'{;
		forval x=1/`nbrands1'{;
			mat aidomega`i'[`x',`x']=(-1+(1+lp1_1[`z',1])*w_nest[`i',`x']
			+lexp1[`z',1]*(1+indexall[`z',1])*w[`i',`x']
			+lpreg_`x'_`x'[`z',1]/bsbar`i'[`x',`z']
			+(laidpreg`x'[`z',1]/bsbar`i'[`x',`z'])*(lp1_1[`z',1]*w_nest[`i',`x']+w[`i',`x']*(1+indexall[`z',1])*lexp1[`z',1]))*x_reg`i'[1,`z']*bsbar`i'[`x',`z'];
		};
	};

	
	
	forval i=1/`nregions'{;
		forval x=1/`nbrands2'{;
			mat aidomega`i'[5+`x',5+`x']=(-1+(1+lp2_2[`z',1])*w_nest[`i',5+`x']
			+lexp2[`z',1]*(1+indexall[`z',1])*w[`i',5+`x']
			+lplight_`x'_`x'[`z',1]/bsbar`i'[5+`x',`z']
			+(laidplight`x'[`z',1]/bsbar`i'[5+`x',`z'])*(lp2_2[`z',1]*w_nest[`i',5+`x']+w[`i',5+`x']*(1+indexall[`z',1])*lexp2[`z',1]))*x_light`i'[1,`z']*bsbar`i'[`x'+5,`z'];
		};
	};

	
	forval i=1/`nregions'{;
		mat aidomega`i'[6,1]=((1+laidpreg1[`z',1]/bsbar`i'[1,`z'])*(lexp1[`z',1]*(1+indexall[`z',1])*w[`i',6]+lp1_2[`z',1]*w_nest[`i',6]))*x_reg`i'[1,`z']*bsbar`i'[1,`z'];
	};
	
	
	
	forval i=1/`nregions'{;
		mat aidomega`i'[1,6]=((1+laidplight1[`z',1]/bsbar`i'[6,`z'])*(lexp2[`z',1]*(1+indexall[`z',1])*w[`i',1]+lp2_1[`z',1]*w_nest[`i',1]))*x_light`i'[1,`z']*bsbar`i'[6,`z'];
	};
	
	
	
	
	forval i=1/`nregions'{;
		mat aidomega`i'[7,2]=((1+laidpreg2[`z',1]/bsbar`i'[2,`z'])*(lexp1[`z',1]*(1+indexall[`z',1])*w[`i',7]+lp1_2[`z',1]*w_nest[`i',7]))*x_reg`i'[1,`z']*bsbar`i'[2,`z'];
	};
	
	
	forval i=1/`nregions'{;
		mat aidomega`i'[2,7]=((1+laidplight2[`z',1]/bsbar`i'[7,`z'])*(lexp2[`z',1]*(1+indexall[`z',1])*w[`i',2]+lp2_1[`z',1]*w_nest[`i',2]))*x_light`i'[1,`z']*bsbar`i'[7,`z'];
	};
	
	
	forval i=1/`nregions'{;
		mat aidomega`i'[8,3]=((1+laidpreg3[`z',1]/bsbar`i'[3,`z'])*(lexp1[`z',1]*(1+indexall[`z',1])*w[`i',8]+lp1_2[`z',1]*w_nest[`i',8]))*x_reg`i'[1,`z']*bsbar`i'[3,`z'];
	};
	
	
	forval i=1/`nregions'{;
		mat aidomega`i'[3,8]=((1+laidplight3[`z',1]/bsbar`i'[8,`z'])*(lexp2[`z',1]*(1+indexall[`z',1])*w[`i',3]+lp2_1[`z',1]*w_nest[`i',3]))*x_light`i'[1,`z']*bsbar`i'[8,`z'];
	};
	
	forval i=1/`nregions'{;
		mat aidomega`i'[9,4]=((1+laidpreg4[`z',1]/bsbar`i'[4,`z'])*(lexp1[`z',1]*(1+indexall[`z',1])*w[`i',9]+lp1_2[`z',1]*w_nest[`i',9]))*x_reg`i'[1,`z']*bsbar`i'[4,`z'];
	};
	
	
	forval i=1/`nregions'{;
		mat aidomega`i'[4,9]=((1+laidplight4[`z',1]/bsbar`i'[9,`z'])*(lexp2[`z',1]*(1+indexall[`z',1])*w[`i',4]+lp2_1[`z',1]*w_nest[`i',4]))*x_light`i'[1,`z']*bsbar`i'[9,`z'];
	};
	
	
	forval i=1/`nregions'{;
		forval x=1/`nbrands1'{;
			mat X`i'[`x',`x']=x_reg`i'[1,`z'];
		};
	};

	forval i=1/`nregions'{;
		forval x=1/`nbrands2'{;
			mat X`i'[`x'+5,`x'+5]=x_light`i'[1,`z'];
		};
	};

	forval i=1/`nregions'{;
		mat aidmu`i'=J(`nbrands',1,0);
		mat aidmu`i'=-inv(aidomega`i')*X`i'*bsbar`i'[1..`nbrands',`z'];		
	};



/*mark-up matrix has 10 columns, which correspond to regional markets, and 7 rows corresponding to brands.  
aidmu are p-mc's with aid system, lmu are p=mc's with linear system, llmu are lerner indices with log-linear system
this just puts them in matrix form as described above*/
	mat aidmu=[aidmu1];
	mat pbar=[bpbar1[1..`nbrands',`z']];
    mat sbar=[bsbar1[1..`nbrands',`z']];
	forval x=2/`nregions'{;
		mat aidmu=[aidmu,aidmu`x'];
		mat pbar=[pbar, bpbar`x'[1..`nbrands',`z']];
		mat sbar=[sbar, bsbar`x'[1..`nbrands',`z']];
	};

	mat pbar=pbar';
	/*mat lpbar=lpbar';*/
	mat sbar=sbar';
	mat aidmu=aidmu';

/*get costs by solving mark-ups.  because element-by-element operations are a pain in stata, i do this in mata*/
/*input markup matrices, average prices and shares, and post-merger ownership matrices into stata*/
	mata: aidmu=st_matrix("aidmu");
	mata: pbar=st_matrix("pbar"); 
	mata: sbar=st_matrix("sbar");
	mata: aidcost=pbar:*(J(`nregions',`nbrands',1)-aidmu);
	mata: st_matrix("aidcost", aidcost);
	

	
    /*defines post-merger ownership matrix of elasticities*/
    mat taidomega=J(`nbrands',`nbrands',0);
	forval i=1/`nregions'{;
		mat taidomega`i'[8,4]=((1+laidpreg4[`z',1]/bsbar`i'[4,`z'])*(lexp1[`z',1]*(1+indexall[`z',1])*w[`i',8]+lp1_2[`z',1]*w_nest[`i',8]))*x_reg`i'[1,`z']*bsbar`i'[4,`z'];
	};

	forval i=1/`nregions'{;
		mat taidomega`i'[4,8]=((1+laidplight3[`z',1]/bsbar`i'[8,`z'])*(lexp2[`z',1]*(1+indexall[`z',1])*w[`i',4]+lp2_1[`z',1]*w_nest[`i',4]))*x_light`i'[1,`z']*bsbar`i'[8,`z'];
	};

	forval i=1/`nregions'{;
		mat taidomega`i'[9,3]=((1+laidpreg3[`z',1]/bsbar`i'[3,`z'])*(lexp1[`z',1]*(1+indexall[`z',1])*w[`i',9]+lp1_2[`z',1]*w_nest[`i',9]))*x_reg`i'[1,`z']*bsbar`i'[3,`z'];
	};


	forval i=1/`nregions'{;
		mat taidomega`i'[3,9]=((1+laidplight4[`z',1]/bsbar`i'[9,`z'])*(lexp2[`z',1]*(1+indexall[`z',1])*w[`i',3]+lp2_1[`z',1]*w_nest[`i',3]))*x_light`i'[1,`z']*bsbar`i'[9,`z'];
	};



	forval i=1/`nregions'{;            
			mat taidomega`i'[4,3]=((1+lp1_1[`z',1])*w_nest[`i',4]
			+lexp1[`z',1]*(1+indexall[`z',1])*w[`i',4]
			+lpreg_4_3[`z',1]/bsbar`i'[3,`z']
			+(laidpreg3[`z',1]/bsbar`i'[3,`z'])*(lp1_1[`z',1]*w_nest[`i',4]+w[`i',4]*(1+indexall[`z',1])*lexp1[`z',1]))*x_reg`i'[1,`z']*bsbar`i'[3,`z'];
	};

	
	forval i=1/`nregions'{;
			mat taidomega`i'[3,4]=((1+lp1_1[`z',1])*w_nest[`i',3]
			+lexp1[`z',1]*(1+indexall[`z',1])*w[`i',3]
			+lpreg_3_4[`z',1]/bsbar`i'[4,`z']
			+(laidpreg4[`z',1]/bsbar`i'[4,`z'])*(lp1_1[`z',1]*w_nest[`i',3]+w[`i',3]*(1+indexall[`z',1])*lexp1[`z',1]))*x_reg`i'[1,`z']*bsbar`i'[4,`z'];
	};



	forval i=1/`nregions'{;
			mat taidomega`i'[9,8]=((1+lp2_2[`z',1])*w_nest[`i',9]
			+lexp2[`z',1]*(1+indexall[`z',1])*w[`i',9]
			+lplight_4_3[`z',1]/bsbar`i'[8,`z']
			+(laidplight3[`z',1]/bsbar`i'[8,`z'])*(lp2_2[`z',1]*w_nest[`i',9]+w[`i',9]*(1+indexall[`z',1])*lexp2[`z',1]))*x_light`i'[1,`z']*bsbar`i'[8,`z'];
	};



	forval i=1/`nregions'{;
		mat taidomega`i'[8,9]=((1+lp2_2[`z',1])*w_nest[`i',8]
			+lexp2[`z',1]*(1+indexall[`z',1])*w[`i',8]
			+lplight_3_4[`z',1]/bsbar`i'[9,`z']
			+(laidplight4[`z',1]/bsbar`i'[9,`z'])*(lp2_2[`z',1]*w_nest[`i',8]+w[`i',8]*(1+indexall[`z',1])*lexp2[`z',1]))*x_light`i'[1,`z']*bsbar`i'[9,`z'];
	};





	#delimit;
	/*solves for post-merger mark-ups*/
	forval i=1/`nregions'{;
		mat postaidomega`i'=aidomega`i'+taidomega`i';
	};


    

	/*solves for post-merger mark-ups*/
    forval i=1/`nregions'{;
        mat define paidmu`i'=J(`nbrands',1,0);
        mat paidmu`i'=-inv(postaidomega`i')*X`i'*bsbar`i'[1..`nbrands',`z'];
    };
    /*puts post-merger mark-ups into a matrix with rows corresponding to markets and columns to products*/
    mat paidmu=[paidmu1];
    forval x=2/`nregions'{;
        mat paidmu=[paidmu,paidmu`x'];
    };
    mat paidmu=paidmu';

    /*now go back to mata to solve for post-merger prices.  will then calculate percentage increases over
    pre-merger average prices for each region, and then take the average over regions for to yield percentage
    price increases for each product.  these average prices are the mxprice variables and the average costs are the
    mxcost variables*/

    mata: paidmu=st_matrix("paidmu");
    mata: invaidp=(J(47,9,1)-paidmu):/aidcost;
    mata: aidp=((J(47,9,1):/invaidp)-pbar):/pbar;
    mata: maidcost=mean(aidcost,1);
    mata: maidprice=mean(aidp,1);
    mata: bsaidcost[`z',.]=maidcost;
    mata: bsaidp[`z',.]=maidprice;
    mata: st_matrix("bsaidcost",bsaidcost);
    mata: st_matrix("bsaidp", bsaidp);
};






/*save nowait*/
#delimit;
drop _all;
svmat bsaidcost;
svmat bsaidp;
save sampledists, replace;

su bsaidcost*, detail;
su bsaidp3, detail;
su bsaidp4, detail;
su bsaidp8, detail;
su bsaidp9, detail;

#delimit;
/*twoway hist bsaidp3 if bsaidp4<1 & bsaidp4>0, title("AID") subtitle("Pennzoil") xtitle("% Price Change in Hundredths") saving(pennapprox, replace);
twoway hist bsaidp4 if bsaidp4<1 & bsaidp4>0, title("AID") subtitle("QuakerState") xtitle("% Price Change in Hundredths") saving(quakerapprox, replace);
*/

#delimit cr;
set more off
mata:
    bp1=st_matrix("bpbar1")
    bp2=st_matrix("bpbar2")
    bp3=st_matrix("bpbar3")
    bp4=st_matrix("bpbar4")
    bp5=st_matrix("bpbar5")
    bp6=st_matrix("bpbar6")
    bp7=st_matrix("bpbar7")
    bp8=st_matrix("bpbar8")
    bp9=st_matrix("bpbar9")
    bp10=st_matrix("bpbar10")
    bp11=st_matrix("bpbar11")
    bp12=st_matrix("bpbar12")
    bp13=st_matrix("bpbar13")
    bp14=st_matrix("bpbar14")
    bp15=st_matrix("bpbar15")
    bp16=st_matrix("bpbar16")
    bp17=st_matrix("bpbar17")
    bp18=st_matrix("bpbar18")
    bp19=st_matrix("bpbar19")
    bp20=st_matrix("bpbar20")
    bp21=st_matrix("bpbar21")
    bp22=st_matrix("bpbar22")
    bp23=st_matrix("bpbar23")
    bp24=st_matrix("bpbar24")
    bp25=st_matrix("bpbar25")
    bp26=st_matrix("bpbar26")
    bp27=st_matrix("bpbar27")
    bp28=st_matrix("bpbar28")
    bp29=st_matrix("bpbar29")
    bp30=st_matrix("bpbar30")
    bp31=st_matrix("bpbar31")
    bp32=st_matrix("bpbar32")
    bp33=st_matrix("bpbar33")
    bp34=st_matrix("bpbar34")
    bp35=st_matrix("bpbar35")
    bp36=st_matrix("bpbar36")
    bp37=st_matrix("bpbar37")
    bp38=st_matrix("bpbar38")
    bp39=st_matrix("bpbar39")
    bp40=st_matrix("bpbar40")
    bp41=st_matrix("bpbar41")
    bp42=st_matrix("bpbar42")
    bp43=st_matrix("bpbar43")
    bp44=st_matrix("bpbar44")
    bp45=st_matrix("bpbar45")
    bp46=st_matrix("bpbar46")
    bp47=st_matrix("bpbar47")

    ppoint=J(1,47,NULL)
    ppoint[1]=&bp1
    ppoint[2]=&bp2
    ppoint[3]=&bp3
    ppoint[4]=&bp4
    ppoint[5]=&bp5
    ppoint[6]=&bp6
    ppoint[7]=&bp7
    ppoint[8]=&bp8
    ppoint[9]=&bp9
    ppoint[10]=&bp10
    ppoint[11]=&bp11
    ppoint[12]=&bp12
    ppoint[13]=&bp13
    ppoint[14]=&bp14
    ppoint[15]=&bp15
    ppoint[16]=&bp16
    ppoint[17]=&bp17
    ppoint[18]=&bp18
    ppoint[19]=&bp19
    ppoint[20]=&bp20
    ppoint[21]=&bp21
    ppoint[22]=&bp22
    ppoint[23]=&bp23
    ppoint[24]=&bp24
    ppoint[25]=&bp25
    ppoint[26]=&bp26
    ppoint[27]=&bp27
    ppoint[28]=&bp28
    ppoint[29]=&bp29
    ppoint[30]=&bp30
    ppoint[31]=&bp31
    ppoint[32]=&bp32
    ppoint[33]=&bp33
    ppoint[34]=&bp34
    ppoint[35]=&bp35
    ppoint[36]=&bp36
    ppoint[37]=&bp37
    ppoint[38]=&bp38
    ppoint[39]=&bp39
    ppoint[40]=&bp40
    ppoint[41]=&bp41
    ppoint[42]=&bp42
    ppoint[43]=&bp43
    ppoint[44]=&bp44
    ppoint[45]=&bp45
    ppoint[46]=&bp46
    ppoint[47]=&bp47
    

        
    bs1=st_matrix("bsbar1")
    bs2=st_matrix("bsbar2")
    bs3=st_matrix("bsbar3")
    bs4=st_matrix("bsbar4")
    bs5=st_matrix("bsbar5")
    bs6=st_matrix("bsbar6")
    bs7=st_matrix("bsbar7")
    bs8=st_matrix("bsbar8")
    bs9=st_matrix("bsbar9")
    bs10=st_matrix("bsbar10")
    bs11=st_matrix("bsbar11")
    bs12=st_matrix("bsbar12")
    bs13=st_matrix("bsbar13")
    bs14=st_matrix("bsbar14")
    bs15=st_matrix("bsbar15")
    bs16=st_matrix("bsbar16")
    bs17=st_matrix("bsbar17")
    bs18=st_matrix("bsbar18")
    bs19=st_matrix("bsbar19")
    bs20=st_matrix("bsbar20")
    bs21=st_matrix("bsbar21")
    bs22=st_matrix("bsbar22")
    bs23=st_matrix("bsbar23")
    bs24=st_matrix("bsbar24")
    bs25=st_matrix("bsbar25")
    bs26=st_matrix("bsbar26")
    bs27=st_matrix("bsbar27")
    bs28=st_matrix("bsbar28")
    bs29=st_matrix("bsbar29")
    bs30=st_matrix("bsbar30")
    bs31=st_matrix("bsbar31")
    bs32=st_matrix("bsbar32")
    bs33=st_matrix("bsbar33")
    bs34=st_matrix("bsbar34")
    bs35=st_matrix("bsbar35")
    bs36=st_matrix("bsbar36")
    bs37=st_matrix("bsbar37")
    bs38=st_matrix("bsbar38")
    bs39=st_matrix("bsbar39")
    bs40=st_matrix("bsbar40")
    bs41=st_matrix("bsbar41")
    bs42=st_matrix("bsbar42")
    bs43=st_matrix("bsbar43")
    bs44=st_matrix("bsbar44")
    bs45=st_matrix("bsbar45")
    bs46=st_matrix("bsbar46")
    bs47=st_matrix("bsbar47")
    

    spoint=J(1,47,NULL)
    spoint[1]=&bs1
    spoint[2]=&bs2
    spoint[3]=&bs3
    spoint[4]=&bs4
    spoint[5]=&bs5
    spoint[6]=&bs6
    spoint[7]=&bs7
    spoint[8]=&bs8
    spoint[9]=&bs9
    spoint[10]=&bs10
    spoint[11]=&bs11
    spoint[12]=&bs12
    spoint[13]=&bs13
    spoint[14]=&bs14
    spoint[15]=&bs15
    spoint[16]=&bs16
    spoint[17]=&bs17
    spoint[18]=&bs18
    spoint[19]=&bs19
    spoint[20]=&bs20
    spoint[21]=&bs21
    spoint[22]=&bs22
    spoint[23]=&bs23
    spoint[24]=&bs24
    spoint[25]=&bs25
    spoint[26]=&bs26
    spoint[27]=&bs27
    spoint[28]=&bs28
    spoint[29]=&bs29
    spoint[30]=&bs30
    spoint[31]=&bs31
    spoint[32]=&bs32
    spoint[33]=&bs33
    spoint[34]=&bs34
    spoint[35]=&bs35
    spoint[36]=&bs36
    spoint[37]=&bs37
    spoint[38]=&bs38
    spoint[39]=&bs39
    spoint[40]=&bs40
    spoint[41]=&bs41
    spoint[42]=&bs42
    spoint[43]=&bs43
    spoint[44]=&bs44
    spoint[45]=&bs45
    spoint[46]=&bs46
    spoint[47]=&bs47
    

/*
    bv1=st_matrix("bvbar1")
    bv2=st_matrix("bvbar2")
    bv3=st_matrix("bvbar3")
    bv4=st_matrix("bvbar4")
    bv5=st_matrix("bvbar5")
    bv6=st_matrix("bvbar6")
    bv7=st_matrix("bvbar7")
    bv8=st_matrix("bvbar8")
    bv9=st_matrix("bvbar9")
    bv10=st_matrix("bvbar10")
    bv11=st_matrix("bvbar11")
    bv12=st_matrix("bvbar12")
    bv13=st_matrix("bvbar13")
    bv14=st_matrix("bvbar14")
    bv15=st_matrix("bvbar15")
    bv16=st_matrix("bvbar16")
    bv17=st_matrix("bvbar17")
    bv18=st_matrix("bvbar18")
    bv19=st_matrix("bvbar19")
    bv20=st_matrix("bvbar20")
    bv21=st_matrix("bvbar21")
    bv22=st_matrix("bvbar22")
    bv23=st_matrix("bvbar23")
    bv24=st_matrix("bvbar24")
    bv25=st_matrix("bvbar25")
    bv26=st_matrix("bvbar26")
    bv27=st_matrix("bvbar27")
    bv28=st_matrix("bvbar28")
    bv29=st_matrix("bvbar29")
    bv30=st_matrix("bvbar30")
    bv31=st_matrix("bvbar31")
    bv32=st_matrix("bvbar32")
    bv33=st_matrix("bvbar33")
    bv34=st_matrix("bvbar34")
    bv35=st_matrix("bvbar35")
    bv36=st_matrix("bvbar36")
    bv37=st_matrix("bvbar37")
    bv38=st_matrix("bvbar38")
    bv39=st_matrix("bvbar39")
    bv40=st_matrix("bvbar40")

    vpoint=J(1,40,NULL)
    vpoint[1]=&bv1
    vpoint[2]=&bv2
    vpoint[3]=&bv3
    vpoint[4]=&bv4
    vpoint[5]=&bv5
    vpoint[6]=&bv6
    vpoint[7]=&bv7
    vpoint[8]=&bv8
    vpoint[9]=&bv9
    vpoint[10]=&bv10
    vpoint[11]=&bv11
    vpoint[12]=&bv12
    vpoint[13]=&bv13
    vpoint[14]=&bv14
    vpoint[15]=&bv15
    vpoint[16]=&bv16
    vpoint[17]=&bv17
    vpoint[18]=&bv18
    vpoint[19]=&bv19
    vpoint[20]=&bv20
    vpoint[21]=&bv21
    vpoint[22]=&bv22
    vpoint[23]=&bv23
    vpoint[24]=&bv24
    vpoint[25]=&bv25
    vpoint[26]=&bv26
    vpoint[27]=&bv27
    vpoint[28]=&bv28
    vpoint[29]=&bv29
    vpoint[30]=&bv30
    vpoint[31]=&bv31
    vpoint[32]=&bv32
    vpoint[33]=&bv33
    vpoint[34]=&bv34
    vpoint[35]=&bv35
    vpoint[36]=&bv36
    vpoint[37]=&bv37
    vpoint[38]=&bv38
    vpoint[39]=&bv39
    vpoint[40]=&bv40
*/  

    laidpreg1=st_matrix("laidpreg1")
    laidpreg2=st_matrix("laidpreg2")
    laidpreg3=st_matrix("laidpreg3")
    laidpreg4=st_matrix("laidpreg4")
    laidpreg5=st_matrix("laidpreg5")
	laidplight1=st_matrix("laidplight1")
    laidplight2=st_matrix("laidplight2")
    laidplight3=st_matrix("laidplight3")
    laidplight4=st_matrix("laidplight4")

    lpreg_1_1=st_matrix("lpreg_1_1")
    lpreg_2_1=st_matrix("lpreg_2_1")
    lpreg_3_1=st_matrix("lpreg_3_1")
    lpreg_4_1=st_matrix("lpreg_4_1")
    lpreg_5_1=st_matrix("lpreg_5_1")
	
	lpreg_1_2=st_matrix("lpreg_1_2")
    lpreg_2_2=st_matrix("lpreg_2_2")
    lpreg_3_2=st_matrix("lpreg_3_2")
    lpreg_4_2=st_matrix("lpreg_4_2")
    lpreg_5_2=st_matrix("lpreg_5_2")
	
	lpreg_1_3=st_matrix("lpreg_1_3")
    lpreg_2_3=st_matrix("lpreg_2_3")
    lpreg_3_3=st_matrix("lpreg_3_3")
    lpreg_4_3=st_matrix("lpreg_4_3")
    lpreg_5_3=st_matrix("lpreg_5_3")
	
	lpreg_1_4=st_matrix("lpreg_1_4")
    lpreg_2_4=st_matrix("lpreg_2_4")
    lpreg_3_4=st_matrix("lpreg_3_4")
    lpreg_4_4=st_matrix("lpreg_4_4")
    lpreg_5_4=st_matrix("lpreg_5_4")
	
	lpreg_1_5=st_matrix("lpreg_1_5")
    lpreg_2_5=st_matrix("lpreg_2_5")
    lpreg_3_5=st_matrix("lpreg_3_5")
    lpreg_4_5=st_matrix("lpreg_4_5")
    lpreg_5_5=st_matrix("lpreg_5_5")
	
	lplight_1_1=st_matrix("lplight_1_1")
    lplight_2_1=st_matrix("lplight_2_1")
    lplight_3_1=st_matrix("lplight_3_1")
    lplight_4_1=st_matrix("lplight_4_1")
	
	lplight_1_2=st_matrix("lplight_1_2")
    lplight_2_2=st_matrix("lplight_2_2")
    lplight_3_2=st_matrix("lplight_3_2")
    lplight_4_2=st_matrix("lplight_4_2")

	lplight_1_3=st_matrix("lplight_1_3")
    lplight_2_3=st_matrix("lplight_2_3")
    lplight_3_3=st_matrix("lplight_3_3")
    lplight_4_3=st_matrix("lplight_4_3")

	lplight_1_4=st_matrix("lplight_1_4")
    lplight_2_4=st_matrix("lplight_2_4")
    lplight_3_4=st_matrix("lplight_3_4")
    lplight_4_4=st_matrix("lplight_4_4")
	
	lp1_1=st_matrix("lp1_1")
	lp1_2=st_matrix("lp1_2")
	lp2_1=st_matrix("lp2_1")
	lp2_2=st_matrix("lp2_2")
	
	lexp1=st_matrix("lexp1")
	lexp2=st_matrix("lexp2")
	indexall=st_matrix("indexall")
	w=st_matrix("w")	
	w_nest=st_matrix("w_nest")	
	
	wnest=st_matrix("wnest")	
	w_all=st_matrix("wall")	
end

drop _all



mata
    e1j=J(1000,9,0)
    e2j=J(1000,9,0)
    e3j=J(1000,9,0)
    e4j=J(1000,9,0)
    e5j=J(1000,9,0)
    e6j=J(1000,9,0)
    e7j=J(1000,9,0)
    e8j=J(1000,9,0)
    e9j=J(1000,9,0)
	

for(z=1;z<=1000;z++){
        /*makes price coefficient matrix for each draw*/
g=(lpreg_1_1[z,1],lpreg_1_2[z,1],lpreg_1_3[z,1],lpreg_1_4[z,1],lpreg_1_5[z,1]
  \lpreg_2_1[z,1],lpreg_2_2[z,1],lpreg_2_3[z,1],lpreg_2_4[z,1],lpreg_2_5[z,1]  
  \lpreg_3_1[z,1],lpreg_3_3[z,1],lpreg_3_3[z,1],lpreg_3_4[z,1],lpreg_3_5[z,1]  
  \lpreg_4_1[z,1],lpreg_4_4[z,1],lpreg_4_3[z,1],lpreg_4_4[z,1],lpreg_4_5[z,1]  
  \lpreg_5_1[z,1],lpreg_5_5[z,1],lpreg_5_3[z,1],lpreg_5_4[z,1],lpreg_5_5[z,1])

g2=(lplight_1_1[z,1],lplight_1_2[z,1],lplight_1_3[z,1],lplight_1_4[z,1]
   \lplight_2_1[z,1],lplight_2_2[z,1],lplight_2_3[z,1],lplight_2_4[z,1]  
   \lplight_3_1[z,1],lplight_3_3[z,1],lplight_3_3[z,1],lplight_3_4[z,1]  
   \lplight_4_1[z,1],lplight_4_4[z,1],lplight_4_3[z,1],lplight_4_4[z,1])
  
  
/*makes expenditure coefficients*/
b=(laidpreg1[z]\laidpreg2[z]\laidpreg3[z]\laidpreg4[z]\laidpreg5[z])
b2=(laidplight1[z]\laidplight2[z]\laidplight3[z]\laidplight4[z])
delta_mid_reg_reg=(lp1_1[z])
delta_mid_reg_light=(lp1_2[z])
delta_mid_light_reg=(lp2_1[z])
delta_mid_light_light=(lp2_2[z])
delta_top=(indexall[z])
beta_reg=(lexp1[z])
beta_light=(lexp2[z])



sx=J(9,47,0)
	for(j=1;j<=47;j++){           
		sx[.,j]=(*spoint[j])[.,z]
	}	            

	s=mean(sx',1)
	s1=s[1..5]'
	s2=s[6..9]'
	
	ereg=J(5,5,0)
	elight=J(4,4,0)
	ereg_light=J(5,4,0)
	elight_reg=J(4,5,0)

for(m=1;m<=5;m++){
	for(n=1;n<=5;n++){
		ereg[m,n]=(1+delta_mid_reg_reg)*wnest[n]+beta_reg*(1+delta_top)*w_all[n]+g[m,n]/s1[m]+(b[m]/s1[m])*(delta_mid_reg_reg*wnest[n]+w_all[n]*(1+delta_top)*beta_reg)
	}
}

for(m=1;m<=4;m++){
	for(n=1;n<=4;n++){
		elight[m,n]=(1+delta_mid_light_light)*wnest[n+5]+beta_light*(1+delta_top)*w_all[n+5]+g2[m,n]/s2[m]+(b2[m]/s2[m])*(delta_mid_light_light*wnest[n+5]+w_all[n+5]*(1+delta_top)*beta_light)
	}
}
ereg=ereg-I(5)
elight=elight-I(4)	



for(m=1;m<=5;m++){
	for(n=1;n<=4;n++){
		ereg_light[m,n]=(1+b[m]/s1[m])*(beta_reg*(1+delta_top)*w_all[n+5]+delta_mid_reg_light*wnest[n+5])
	}
}

for(m=1;m<=4;m++){
	for(n=1;n<=5;n++){
		elight_reg[m,n]=(1+b2[m]/s2[m])*(beta_light*(1+delta_top)*w_all[n]+delta_mid_light_reg*wnest[n])
	}
}


		e1j[z,1]=ereg[1,1]
		e1j[z,2]=ereg[1,2]
		e1j[z,3]=ereg[1,3]
		e1j[z,4]=ereg[1,4]
		e1j[z,5]=ereg[1,5]
		e1j[z,6]=ereg_light[1,1]
		e1j[z,7]=ereg_light[1,2]
		e1j[z,8]=ereg_light[1,3]
		e1j[z,9]=ereg_light[1,4]
		
		e2j[z,1]=ereg[2,1]
		e2j[z,2]=ereg[2,2]
		e2j[z,3]=ereg[2,3]
		e2j[z,4]=ereg[2,4]
		e2j[z,5]=ereg[2,5]
		e2j[z,6]=ereg_light[2,1]
		e2j[z,7]=ereg_light[2,2]
		e2j[z,8]=ereg_light[2,3]
		e2j[z,9]=ereg_light[2,4]
		
		e3j[z,1]=ereg[3,1]
		e3j[z,2]=ereg[3,2]
		e3j[z,3]=ereg[3,3]
		e3j[z,4]=ereg[3,4]
		e3j[z,5]=ereg[3,5]
		e3j[z,6]=ereg_light[3,1]
		e3j[z,7]=ereg_light[3,2]
		e3j[z,8]=ereg_light[3,3]
		e3j[z,9]=ereg_light[3,4]
		
		e4j[z,1]=ereg[4,1]
		e4j[z,2]=ereg[4,2]
		e4j[z,3]=ereg[4,3]
		e4j[z,4]=ereg[4,4]
		e4j[z,5]=ereg[4,5]
		e4j[z,6]=ereg_light[4,1]
		e4j[z,7]=ereg_light[4,2]
		e4j[z,8]=ereg_light[4,3]
		e4j[z,9]=ereg_light[4,4]
		
		e5j[z,1]=ereg[5,1]
		e5j[z,2]=ereg[5,2]
		e5j[z,3]=ereg[5,3]
		e5j[z,4]=ereg[5,4]
		e5j[z,5]=ereg[5,5]
		e5j[z,6]=ereg_light[5,1]
		e5j[z,7]=ereg_light[5,2]
		e5j[z,8]=ereg_light[5,3]
		e5j[z,9]=ereg_light[5,4]
		
		e6j[z,1]=elight_reg[1,1]
		e6j[z,2]=elight_reg[1,2]
		e6j[z,3]=elight_reg[1,3]
		e6j[z,4]=elight_reg[1,4]
		e6j[z,5]=elight_reg[1,5]
		e6j[z,6]=elight[1,1]
		e6j[z,7]=elight[1,2]
		e6j[z,8]=elight[1,3]
		e6j[z,9]=elight[1,4]
		
		e7j[z,1]=elight_reg[2,1]
		e7j[z,2]=elight_reg[2,2]
		e7j[z,3]=elight_reg[2,3]
		e7j[z,4]=elight_reg[2,4]
		e7j[z,5]=elight_reg[2,5]
		e7j[z,6]=elight[2,1]
		e7j[z,7]=elight[2,2]
		e7j[z,8]=elight[2,3]
		e7j[z,9]=elight[2,4]
		
		e8j[z,1]=elight_reg[3,1]
		e8j[z,2]=elight_reg[3,2]
		e8j[z,3]=elight_reg[3,3]
		e8j[z,4]=elight_reg[3,4]
		e8j[z,5]=elight_reg[3,5]
		e8j[z,6]=elight[3,1]
		e8j[z,7]=elight[3,2]
		e8j[z,8]=elight[3,3]
		e8j[z,9]=elight[3,4]
		
		e9j[z,1]=elight_reg[4,1]
		e9j[z,2]=elight_reg[4,2]
		e9j[z,3]=elight_reg[4,3]
		e9j[z,4]=elight_reg[4,4]
		e9j[z,5]=elight_reg[4,5]
		e9j[z,6]=elight[4,1]
		e9j[z,7]=elight[4,2]
		e9j[z,8]=elight[4,3]
		e9j[z,9]=elight[4,4]
		
		
}
mean(e1j)
mean(e2j)
mean(e3j)				      
mean(e4j)
mean(e5j)
mean(e6j)
mean(e7j)				      
mean(e8j)
mean(e9j)

aidsolselas=J(9,9,0)

aidsolselas[1,.]=sqrt(mean(e1j:*e1j)-mean(e1j):*mean(e1j))				      
aidsolselas[2,.]=sqrt(mean(e2j:*e2j)-mean(e2j):*mean(e2j))				      
aidsolselas[3,.]=sqrt(mean(e3j:*e3j)-mean(e3j):*mean(e3j))
aidsolselas[4,.]=sqrt(mean(e4j:*e4j)-mean(e4j):*mean(e4j))
aidsolselas[5,.]=sqrt(mean(e5j:*e5j)-mean(e5j):*mean(e5j))				      
aidsolselas[6,.]=sqrt(mean(e6j:*e6j)-mean(e6j):*mean(e6j))				      
aidsolselas[7,.]=sqrt(mean(e7j:*e7j)-mean(e7j):*mean(e7j))				      
aidsolselas[8,.]=sqrt(mean(e8j:*e8j)-mean(e8j):*mean(e8j))				      
aidsolselas[9,.]=sqrt(mean(e9j:*e9j)-mean(e9j):*mean(e9j))				      				      

end





cap log close

log using syrelasticities, replace
	/*table with caption "Syrup Elasticities, AIDS Model with OLS" has 10 rows and 5 columns. The matrix below contains standard errors for 
	these elasticities.  The first row of this matrix goes in the second row of "Syrup Elasticities, AIDS Model with OLS", the second row of this matrix 
	goes in the 4th, etc."*/
	mata: aidsolselas
log close



set more off
mata
    bsfaid=J(1000,9,0)
	bsfaidfocs=J(1000,9,0)
	regsolution=J(1000,1,0)
    for(z=1;z<=1000;z++){
		g=(lpreg_1_1[z,1],lpreg_1_2[z,1],lpreg_1_3[z,1],lpreg_1_4[z,1],lpreg_1_5[z,1]
		\lpreg_2_1[z,1],lpreg_2_2[z,1],lpreg_2_3[z,1],lpreg_2_4[z,1],lpreg_2_5[z,1]  
		\lpreg_3_1[z,1],lpreg_3_3[z,1],lpreg_3_3[z,1],lpreg_3_4[z,1],lpreg_3_5[z,1]  
		\lpreg_4_1[z,1],lpreg_4_4[z,1],lpreg_4_3[z,1],lpreg_4_4[z,1],lpreg_4_5[z,1]  
		\lpreg_5_1[z,1],lpreg_5_5[z,1],lpreg_5_3[z,1],lpreg_5_4[z,1],lpreg_5_5[z,1])

		g2=(lplight_1_1[z,1],lplight_1_2[z,1],lplight_1_3[z,1],lplight_1_4[z,1]
		\lplight_2_1[z,1],lplight_2_2[z,1],lplight_2_3[z,1],lplight_2_4[z,1]  
		\lplight_3_1[z,1],lplight_3_3[z,1],lplight_3_3[z,1],lplight_3_4[z,1]  
		\lplight_4_1[z,1],lplight_4_4[z,1],lplight_4_3[z,1],lplight_4_4[z,1])
  
  
		/*makes expenditure coefficients*/
		b=(laidpreg1[z]\laidpreg2[z]\laidpreg3[z]\laidpreg4[z]\laidpreg5[z])
		b2=(laidplight1[z]\laidplight2[z]\laidplight3[z]\laidplight4[z])
		delta_mid_reg_reg=(lp1_1[z])
		delta_mid_reg_light=(lp1_2[z])
		delta_mid_light_reg=(lp2_1[z])
		delta_mid_light_light=(lp2_2[z])
		delta_top=(indexall[z])
		beta_reg=(lexp1[z])
		beta_light=(lexp2[z])
		delta=(delta_mid_reg_reg, delta_mid_reg_light\delta_mid_light_reg, delta_mid_light_light)
		beta=(beta_reg\beta_light)

		temp=J(47,9,0)
		check=J(47,1,0)	
		aidfocs=J(47,9,0)
        for(j=1;j<=47;j++){
            p=(*ppoint[j])[.,z]
            pinit=(*ppoint[j])[.,z]
            s=(*spoint[j])[.,z]
            aidomega=J(9,9,0)
		/*******************************START HERE*******************************/
			for(x=1;x<=5;x++){
				aidomega[x,x]=(-1+(1+delta_mid_reg_reg)*w_nest[j,x]+beta_reg*(1+delta_top)*w[j,x]+g[x,x]/s[x]+(b[x]/s[x])*(delta_mid_reg_reg*w_nest[j,x]+w[j,x]*(1+delta_top)*beta_reg))*X[j,1]*s[x]
			}
			for(x=6;x<=9;x++){
				aidomega[x,x]=(-1+(1+delta_mid_light_light)*w_nest[j,x]+beta_light*(1+delta_top)*w[j,x]+g2[x-5,x-5]/s[x]+(b2[x-5]/s[x])*(delta_mid_light_light*w_nest[j,x]+w[j,x]*(1+delta_top)*beta_light))*X[j,2]*s[x]
			}
			aidomega[6,1]=((1+b[1]/s[1])*(beta_reg*(1+delta_top)*w[j,6]+delta_mid_reg_light*w_nest[j,6]))*X[j,1]*s[1]			            			
			aidomega[1,6]=((1+b2[1]/s[6])*(beta_light*(1+delta_top)*w[j,1]+delta_mid_light_reg*w_nest[j,1]))*X[j,2]*s[6]
			aidomega[7,2]=((1+b[2]/s[2])*(beta_reg*(1+delta_top)*w[j,7]+delta_mid_reg_light*w_nest[j,7]))*X[j,1]*s[2]
			aidomega[2,7]=((1+b2[2]/s[7])*(beta_light*(1+delta_top)*w[j,2]+delta_mid_light_reg*w_nest[j,2]))*X[j,2]*s[7]
			aidomega[8,3]=((1+b[3]/s[3])*(beta_reg*(1+delta_top)*w[j,8]+delta_mid_reg_light*w_nest[j,8]))*X[j,1]*s[3]
			aidomega[3,8]=((1+b2[3]/s[8])*(beta_light*(1+delta_top)*w[j,3]+delta_mid_light_reg*w_nest[j,3]))*X[j,2]*s[8]
			aidomega[9,4]=((1+b[4]/s[4])*(beta_reg*(1+delta_top)*w[j,9]+delta_mid_reg_light*w_nest[j,9]))*X[j,1]*s[4]
			aidomega[4,9]=((1+b2[4]/s[9])*(beta_light*(1+delta_top)*w[j,4]+delta_mid_light_reg*w_nest[j,4]))*X[j,2]*s[9]
			
			Xi=J(9,9,0)
			Xi[1,1]=X[j,1]
			Xi[2,2]=X[j,1]
			Xi[3,3]=X[j,1]
			Xi[4,4]=X[j,1]
			Xi[5,5]=X[j,1]
			Xi[6,6]=X[j,2]
			Xi[7,7]=X[j,2]
			Xi[8,8]=X[j,2]
			Xi[9,9]=X[j,2]

			mu=-luinv(aidomega)*Xi*s
			aidcost=p:*(J(9,1,1)-mu)
			pinit=(*ppoint[j])[.,z]
		
			a=csolve(&msyraids(),pinit,1e-6,400,s,p,aidcost,g,g2,b,b2,w[j,.]',w_nest[j,.]',X[j,.]',delta,beta,delta_top)		
			ws=msyraids(a[1..9,.],s,p,aidcost,g,g2,b,b2,w[j,.]',w_nest[j,.]',X[j,.]',delta,beta,delta_top)
			check[j,1]=((a[10,1]+(a[1..rows(a)-1,.]==J(rows(a)-1,1,.))+(ws[1,.]==.)+(ws[2,.]==.)+(ws[3,.]==.)+(ws[4,.]==.)+(ws[5,.]==.)+(ws[6,.]==.)+(ws[7,.]==.)+(ws[8,.]==.)+(ws[9,.]==.))==0)
			a=a[1..9,.]
			a=(a-p):/p
			temp[j,.]=a'
			aidfocs[j,.]=msyraids(postdirect[.,j],s,p,aidcost,g,g2,b,b2,w[j,.]',w_nest[j,.]',X[j,.]',delta,beta,delta_top)'
    }
			bsfaid[z,.]=mm_median(select(temp,check))
			bsfaidfocs[z,.]=mean(select(aidfocs,check))
}       
	st_matrix("bsfaid",bsfaid)
	st_matrix("bsfaidfocs", bsfaidfocs)
end


svmat bsfaidfocs
svmat bsfaid

log using syrtable2col3, replace
centile bsfaid1, centile(2.5, 97.5)
centile bsfaid2, centile(2.5, 97.5)
centile bsfaid3, centile(2.5, 97.5)
centile bsfaid4, centile(2.5, 97.5)
centile bsfaid5, centile(2.5, 97.5)
centile bsfaid6, centile(2.5, 97.5)
centile bsfaid7, centile(2.5, 97.5)
centile bsfaid8, centile(2.5, 97.5)
centile bsfaid9, centile(2.5, 97.5)

centile bsfaidfocs1, centile(2.5, 97.5)
centile bsfaidfocs2, centile(2.5, 97.5)
centile bsfaidfocs3, centile(2.5, 97.5)
centile bsfaidfocs4, centile(2.5, 97.5)
centile bsfaidfocs5, centile(2.5, 97.5)
centile bsfaidfocs6, centile(2.5, 97.5)
centile bsfaidfocs7, centile(2.5, 97.5)
centile bsfaidfocs8, centile(2.5, 97.5)
centile bsfaidfocs9, centile(2.5, 97.5)

log close
save c:/data/restat_12_9/bootstrap_results/syraidsols, replace

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

forval x=1/9{
	replace bsfaid`x'=100*bsfaid`x'
}

rename bsfaid1 Aunt_Jemima
rename bsfaid2 Hungry_Jack
rename bsfaid3 Log_Cabin
rename bsfaid4 Mrs_Butterworth
rename bsfaid5 PrivateLabel

xxxxxxxxxxxxxxxxxxx
twoway histogram Castrol, title("Castrol") xtitle("Percent Price Change") saving(hc, replace)
twoway hist Havoline, title("Havoline") xtitle("Percent Price Change") saving(hh, replace)
twoway histogram Mobil, title("Mobil") xtitle("Percent Price Change") saving(hm, replace)
twoway histogram Pennzoil, title("Pennzoil") xtitle("Percent Price Change") saving(hpe, replace)
twoway hist PrivateLabel, title("Private Label") xtitle("Percent Price Change") saving(hp, replace)
twoway hist Quaker, title("Quaker") xtitle("Percent Price Change") saving(hq, replace)
twoway hist Valvoline, title("Valvoline") xtitle("Percent Price Change") saving(hv, replace)

graph combine hc.gph hh.gph hm.gph hpe.gph hp.gph hq.gph hv.gph, title(Sampling Distributions of Price Changes)
graph export olsaidsoil.eps, replace

/*Are sampling distributions normal?*/
/*qq plots*/
qnorm Castrol, title("Castrol") saving(qc, replace)
qnorm Havoline, title("Havoline") saving(qh, replace)
qnorm Mobil, title("Mobil") saving(qm, replace)
qnorm Pennzoil, title("Pennzoil") saving(qp, replace)
qnorm PrivateLabel, title("Private Label") saving(qplo, replace)
qnorm Quaker, title("Quaker State") saving(qq, replace)
qnorm Valvoline, title("Valvoline") saving(qv, replace)


graph combine qc.gph qh.gph qm.gph qp.gph qplo.gph qq.gph qv.gph
graph export qolsaidsoil.eps, replace

/*Jarque-Bera Tests for Normality*/
egen sc=skew(Castrol)
egen kc=kurt(Castrol)
gen jbc=_N*(saj^2/6+(kaj-3)^2/24)

egen sh=skew(Havoline)
egen kh=kurt(Havoline)
gen jbh=_N*(sh^2/6+(kh-3)^2/24)

egen sm=skew(Mobil)
egen km=kurt(Mobil)
gen jbm=_N*(sm^2/6+(km-3)^2/24)

egen sp=skew(Pennzoil)
egen kp=kurt(Pennzoil)
gen jbp=_N*(sp^2/6+(kp-3)^2/24)

egen splo=skew(PrivateLabel)
egen kplo=kurt(PrivateLabel)
gen jbplo=_N*(splo^2/6+(kplo-3)^2/24)

egen sq=skew(QuakerState)
egen kq=kurt(QuakerState)
gen jbq=_N*(sq^2/6+(kq-3)^2/24)

egen sv=skew(Valvoline)
egen kv=kurt(Valvoline)
gen jbv=_N*(sv^2/6+(kv-3)^2/24)

sum jbaj jbhj jblc jbmb jbplo
save oilaidsols, replace









xxxxx
